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Spectral Measurements of Shock Layer Radiation 
in an Arc-Jet Wind Tunnel 

Giuseppe Palumbo , 1 Roger Craig 2 and Armando Carrasco 2 
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ABSTRACT 

Measurements were made of the radiating gas cap of a blunt body in an NASA- 
Ames 20 MW arc-jet wind tunnel. The test gas was air. Spectra of the flux 
incident on a small aperture centered at the stagnation region were obtained. A 
helium-cooled, magnesium fluoride window transmuted flux into an evacuated 
collimating system that ocused the aperture onto the entrance slit of a 
spectrometer. Data were obtained with films and by photomultipliers. The 
range covered was 1 20 nm to 1 000 nm and the resolution was 0.05 nm to 0.5 nm. 
This paper presents preliminary results from the experiment. Representative 
spectral records from 200 nm to 1 000 nm are shown. The spectra show the 
atomic lines from oxygen and nitrogen in the IR, as well as the molecular 
systems of NO, N2, Nz + , and CN. Copper, as a contaminant, and carbon are 
tentatively identified. Planned subsequent laboratory work will result in 
calibration of spectral sensitivity and refined wavelength determinations. 
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INTRODUCTION 

Arc-jet wind tunnels, simulating conditions of nigh speed flight, are likely to 
be more and more important for aerospace research in the future. 

Computational modelling -eplacing the work of low speed wind tunnels, wherein 
real gas effects do not dominate, is becoming a reality. However, high speed 
flights, such as earth entry from a space mission (low earth orbit, LEO, or 
further), involves fiowfieid gases which are not in thermochemical equilibrium 
(references 1 and 2). This can include important levels of transport of radiative 
energy (reference 3). At present we cannot accurately predict the 
environmental conditions surrounding these spacecraft. Figure I shows flight 
regimes experienced for the entries of the Shuttle, Apollo and a proposed 
Aeroassisted Space Transfer Vehicle (ASTV). An ASTV is a conceptual vehicle 
which utilizes aerodynamic forces to decelerate and alter orbit parameters to 
rendezvous with the Shuttle or a space station (reference 1 ). The design of 
future, highly efficient heat-shields (i. e., not overdesigned) for these vehicles, 
and others, requires the capability of making accurate heating predictions. This 
fact is acknowledged by the inception of the Aeroassist Flight Experiment 
(reference 4). This experiment utilized a vehicle designed to be deployed and 
recovered by the Shuttle. It was to fly a trajectory which simulated an ASTV 
maneuver from geosynchronous earth orbit (GEO) to LEO for rendezvous with the 
Shuttle or a space station. The objective was to gather a data base for 
development of the ASTV and other advanced space transportation systems. The 
program was recently canceled due to funding limitations. Although the 
activity was cancelled the technical requirement remains. 

Real gas computational models are being developed (reference 5) but validating 
data are lacking. Various models predict widely different radiative heating 
levels For example, predictions of ASTV radiative heating from only the VUV 
atomic lines range from insignificant amounts to levels dominating the overall 
heating (references 6, 7 and 8). 

Arc-jet wind tunnels can produce the enthalpy and pressure conditions 
simulating these high speed entries. These facilities generally consist of an 
arc heater, a supersonic nozzle, a test box, a model holder and the necessary 
equipment to deal with the extreme heat transfer and the exhaust Reference 9 
describes arc-jet wind tunnels in some detail. In the past their main use has 
been for heat shield development. Future uses will include aerothermodynamic 
testing for spacecraft flight. Also these facilities can be used to conduct 
experiments to validate computational codes. Even though the size scale cannot 
be simulated, a flow-field rich in non-equiiibrated, radiating gas can be 


generated. Computational models can be exercised on arc -jet test conditions 
for comparison with experimental results. There are some major difficulties in 
this approach, however. The free stream plasma flow conditions (enthalpy, 
species distribution, energy states distribution, etc.) in an arc -jet are not well 
understood. Initial conditions for shock layer computational models can only be 
estimated. At present there are efforts to calculate the moael test 
environments using computational models starting with the arc column, 
continuing with the conical or contoured expansion nozzle, and culminating in 
the flow field of a model located in the exit flow from the nozzle (reference 
10). Experimental results are needed to support the theoretical work. A review 
of techniques used to study arc-jet wind tunnel flow is given in reference 1 1 . 

This paper reports on initial results of an arc -jet wind tunnel experiment. The 
experiment measured the spectral radiative flux emanating from the shock la ye r 
and incident on a blunt model placed in the supersonic stream. Data was taken 
from 1 20 nm to 900. This paper presents resuits the longer wavelength resuits 
from 200 nm. The VUV results will be reported on elsewhere. The purpose was 
to gather data to help characteriz - the stagnation region shock layer flow. 
Understanding this flow would be valuable for development and validation of 
advanced arc-jet wind tunnel design codes and help extend the capability of this 
class of facility to advanced aerothermodynamic testing. 

The data shown in this paper are preliminary and do not benefit from final 
calibrations. The flux was spectrally resolved and over a wide range. The data 
will be analyzed and used to identify important radiating species and help 
determine the state of the gas and the operating characteristics of arc-jet wind 
tunnels. 


EXPERIMENTAL SETUP 

The experimental setup is as described in reference 1 2 except that the model 
diameter was increased to 6 inches. Figures 2 and 3 are from this reference and 
are included here. Figure 2 is a schematic plan view of the experiment. 
Supersonic flow from the arc is produced in the nozzle and the model is placed 
in the free stream. The standing shock layer over the model is indicated. Figure 
3 is a schematic view of the model. The model face is a 6" diameter flat disk 
and the entire model is water cooled. A small aperture centered on the forward 
face admits *he surface radiative flux. A helium cooled window, shown 
immersed in a cavity below the apei cure, transmits the surface flux into the 
optical system. The optical system consists of a flat mirror, directing th 
optical path to a concave mirror, which in turn collimates the beam to another 
concave mirror shown at the top of figure 2. This mirror focuses the aperture, 
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vis a flat turning miiror, onto the slit of the spectrograph. The spectrograph 
was a 0.5 meter instrument using a modified Czerny-Turner optical layout with 
the provision for operating either as a scanning monochromator or a fiim 
spectrograph . The portions of the system exposed to the plasma were ali water 
cooied. Figure 4 shows the model arrangement in the test box. The nozzle exit 
can be seen at the left, and the apertu r e can be seen centered on the model face. 
The optical axis in front of the me del was canted 1 5° from the centerline to 
reject radiation from the arc column. 

TEST CONDITIONS 

The test was conducted in the Ames 20 Megawatt Aerodynamic Heating Facility. 
The facility was operated with a supersonic nozzle with a 1.5 inch diameter 
throat and an 1 8" exit diameter. The facility operating parameters were as 
follows: 


Test gas mixture: 80% air" and 20% Argon by mass, 

arc current: 1 000 amperes, 

arc column pressure: 1 5 psia, 

nozzle pressure: 7.5 psia, 

test box pressure: 0.2 to 0.3 mm Hg, and 

stagnation pressure: 9 mm Hg. 

The resulting free stream conditions for the present experiment is approximate 
as indicated in by the circle in figure 1 . The velocity is deduced from an 
estimate of the free stream enthalpy which was, in turn, derived from the 
heating rate of a small sphere placed in the stream, it should be noted that the 
arc-jet free stream enthalpy, although repeatable, is not well characterized. 
Indeed, one of the objectives of this experiment is to help quantify this and 
other arc-jet wind tunnel performance characteristics. 

Figure 5 is a photograph taken during a test. The shock is seen well formed over 
the flat model face. The intense radiation from the high temperature shock 
layer gases is clearly evident. 

The experiment consisted of a series of 1 0 minute runs during each of which a 
specific spectral region was examined. The length of a run was limited by the 
heating loads on the test box. Data were obtained using film and 
photomultiplier tubes. The photomultiplier signal was processed with 
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synchronous amplification. Readout was on a cha recorder and diaitaily 
recorded. 

The instrument configurations used in the tests are summarized in Table 1 . 

DIGITAL DATA ACQUISITION 


Four signals were recorded digitally during the experiment tests: stagnation 
pressure on the mode! face, optical system vacuum pressure, photomultiplier 
high voitage, and photomultiplier output signs!. The stagnation pressure, vacuum 
pressure, and the high voitage were recorded at one hertz with an analog mixer. 
The photomultiplier output signal was recorded ‘crty hertz. 

Figure 6 shows the configuration of the data recording system for the time 
history of the pressure, vacuum, and high voltage signals. The system used a 
multiplexer and a controller cr nfigured with a portable personal computer. 
Twisted and shielded grounded cables were used to transmit the signals. The 
data acquisition and instrument control software were configured to record the 
data during each run at one hertz and concurrently display it on the monitor. 

Figuie 7 shows the configuration of the data recording system for the 
photomultiplier records. The system consisted of a controller, a scanning A/D 
converter, and a 1 6 channel isolated input rack connected to an NB-DMA board 
contained within the personal computer system. The data acquisition and 
control software was set to record at 40 Hertz. 

The photomultiplier data presented here were manipulated for display purposes 
using graphic and analysis software. 

RESULTS AND DISCUSSION 

Computed Spectral Details. 

The development of computational codes to predict the shock layer radiation, 
such as that measured by this experiment, is a parallel, ongoing activity at the 
Ames Research Center. These codes are based on calculations of the arc 
process, the expansion in the nozzle and the shock la'"" processes. The 
radiating shock region is very non-homogeneous and in .olves many kinetically 
controlled processes. Substantial radiation emanates from the nonequilibrated 
regions immediately behind the shock. Here the kinetic temperature is 
extremely high, approaching 50,000K. The spectral line radiation is strongly 
Doppier broadened by this high temperature. The cooler parts of the shock will 


partially absorb the strong features directed toward the surface, but absorption 
will occur only at around the line centers because the absorption lines are 
narrow due to the low kinetic temperature. The surface radiative heating flux 
is thus seen to be the aggregate radiation from regions of very different 
conditions and estimates from computational codes involve consideration of 
complex interactive and kinetic processes. At present these codes exist as 
separate codes but the goal is to couple them together and validate the result as 
a reliable code for predicting arc-jet wind tunnel conditions. Figure 8 is a 
calculated spectrum. NO bands are seen from 1 90 nm to 300 nm. N2(2+) bands 
contribute energy from 280 nm on and merge with the stronger N2 + (1-) and 
N 2 ( 1 -r- ) bands which dominate the molecular systems from 300 nm to the 
infrared. Rich and diagnostically important atomic line radiation is predicted 
from oxygen and nitrogen in the infrared. 

Test Results 

As indicated above the data were taken with turn and with photomultiplier 
tubes. The data presented herein are not corrected for instrument spectral 
response. This correction will be done with future work. The photomultiplier 

data was digitized in-situ. The film spectra density will be subsequently 
digitized. M y 

Ultraviolet-visible results. 


Data was acquired in the spectral region from 200 nm to 500 nm by both film 
and photomultiplier tube. Figure 9 is a print from a film record taken on 
8 1 0/92. The spectral range is from 330 nm to 450 nm. The uppermost 
spectrum is from the shock layer radiation and includes N 2 + (l-) B 2 Z u + -X 2 Z g + 
and CNv B 2 Z~X 2 Z band systems and atomic lines from nitrogen and oxygen The 
other spectra are mercury line calibrations. Figure 1 0 is a densitometer scan of 
the 360 nm to 395 nm region of figure M in which details of a portion of the 

Nz+O-) system, are shown. The vib-rot structure can be seen in detail. This and 
the other film records wiil be digitally analyzed for subsequent use. 

Photomultiplier data were obtained in this region also. Figure 1 1 shows 
photomultiplier data of the spectral details from 200 nm to 300 nm. This 
spectral region is dominated by the NOy A 2 r-X 2 n bands. Figure 1 2 is the 
photomultiplier data from 250 nm to 500 nm. These data were recorded 
digitally in-situ as well as with chart recorder backup. As with the VUV data, 
there is no correction 'o these figures to account for instrument spectral 
response. This correction awaits further laboratory work. 


Vcible-infrared resuits. 


f *’ ' i: 2 limitations prohibited use of photomultiplier tubes beyond 700 nm. Fiim 
s.. ctra, however, were obtained with great detail. Figure 13 is a print of the 
fil l taken from the 7/23/92 test, it covers the spectral range from 450 nm to 
950 nm with a resolution of approximately 0.1 nm. Figure 14 is a densitometer 
trace of this record. The details which overlap with the photon ultiplier record 
can be seen. The important atomic oxygen and nitrogen emission lines are well 
resolved. Figure 1 5 is an densitometer trace of an expanded portion of this 
fiim indicating the resolution of the atomic lines. Several features await 
identification. Calibrations of the non-linearities of the film system will 
result in a data base allowing comparison of the intensities of these features 
with theoretical resuits. Their relative intensities will help understand the 
electron density of the plasma. 


CONCLUSIONS 

The spectra have been obtained from the flux incident on the stagnation surface 
of a fiat mode! placed in the supersonic stream in an arc -jet wind tunnel. The 
preliminary data set shown is detailed at high resolution from 200 nm to 1000 
nm. Important radiators are evident. The data set will be further developed and 
used to he'p refine and calibrate computational models of the 
aerothermodynamics of entry and of arc-jet wind tunnels. 
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TABLE 1 

l-TESTS USING FILM 


Date 

Grating 

Blaze 

Slit 

Exposure, 

X range 


Ruling, 

l/mm 

Wavelength 

nm 

microns, 

seconds 

nm 

7/8/92 

60 

300 

20 

100 

200 to 500 

7/15/92 

180 

150 

20 

60, 300 

200 to 500 

7/1 6/92 

300 

500 

20 

10, 50, 300 

300 to 1000 

7/23/92 

300 

500 

20 

100, 390 

300 to 1000 

7/24/92 

1200 

150 

20 

100, 600 

200 to 400 

8/10/92 

1200 

150 

20 

10, 100, 480 

300 to 500 


II-TESTS USING PHOTOMULTIPLIER TUEE 



Grating 

Blaze 

resolution, 

spectral 

X range, 


l/mm 

nm 

nm 

features 

nm 

8/4/92 

2400 

150 

0.4 

N2+(1-) 

3C0 to 500 

8/4/92 

2400 

150 

0.1 

N2+(1-) 

300 to 500 

8/5/92 

2400 

150 

0.4 

NO/ 

200 to 290 

8/5/92 

2400 

150 

0.1 

NOy 

200 to 290 
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FIGURES 


Figure 1 . 

Might regimes of mission returns for ASTV, Shuttle, and Apollo. The point 
indicated by the circle is an estimate of the altitude and velocity simulated by 
the a c-jet conditions used during the present tests. 

Figure 2. 

Schematic diagram of experimental setup. The arc column is to the left of the 
nozzle. 


Figure 3. 

Schematic diagram of model showing the aperture and window and orientation 
of the bow shock. Radiation incident onto the surface is reflected by the 
turning mirror to the concave mirror and thence to the spectrograph. Vacuum 
was maintained <0.0 Ip for the VUV tests. 

Figure 4. 

Photograph of model in test box. The optical system is protected by the water 
cooled coils shown as well as internal water cooling. The aperture can be seen 
in the face of the model. The nozzle exit can be seen at the left. 

Figure 5. 

Photograph of model during test. The stagnation shock is seen well formed over 
the model face and the intense radiation from the shock heated air is evident. 

Figure 6. 

Data acquisition system block diagram for recording the stagnation pressure on 
the model face, optical system vacuum pressure, photomultiplier high voltage. 

Figure 7. 

Data acquisition system block diagram for recording the photomultiplier output 
signal. 


Figure 8. 

Preliminary calculation of the spectral surface flux from the stagnation region 
of a model placed in an arc-jet wind tunnel. 


Figure 9 

Print fr°m a film record taken on 8/10/92. The spectral range is from 330 nm 
lO 450 nm. , he uppermost spectrum is from the shock layer radiation and 

TTT “/‘Ig and Cfiv S2I-X2I band systems and atomic lines 

" om nitr °geh and oxygen. The other spectra are mercury line calibrations. 

Figure 10 

Densitometer trace of the 360 nm to 395 nm region of the film record of figure 
VI is shown. This region includes a portion of the N?+(l-) B 2 I U ~-X 2 ^V band 
system. The vibration-rotation lines are seen to be'wefl resoled. The atomic 
oxygen line at 394.8 nm is also shown. 


Figure 1 1 

Photomultiplier record of spectrum from 200 nm to 300 nm. This spectral 
region is dominated by the NO-y A 2 r -X 2 TI bands. 

Figure 1 2. 

Photomi/tiplier record of spectrum from 250 nm to 500 nm. The molecular 
spectra from NO-, A2r-X2n, N 2 +(1-) 82Zg*-x2Z,f and CNv B2X-X2Z are in this 
region as well as atomic lines from oxygen and nitrogen and copper 
(contamination from the arc electrodes). 

nu ... Figure 1 3. 

Photographic print of 7/23/92 film spectrum from 450 nm to 950 nm. Two 
different duration exposures of the shock layer emission were taken. Exposure 
A is a short exposure of emission from the shock layer and an overlapping 
calibration exposure from a Hg lamp. Exposures B and D are Hg lamp calibration 
exposures. Exposure C is a long exposure without overlapping Hg lamp 

calibration. Selected Hg emission lines in the first, second and third orders are 
indicated. 


Figure 14. 

Densitometer trace of 7/23/92 film spectrum from 500 nm to 1 000 nm. The 
important atomic oxygen and nitrogen emission lines are well resolved. Also 

seen are the copper lines at 324.7 nm and 327.4 nm in the second and third 
orders. 


Figure 1 5. 

Expanded portion of the 7/23/92 film spectrum densitometer trace showing the 
atomic lines near between 800 nm and 880 nm. Several features await 
identification. 
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ABSTRACT 


The FL3D solver, a full three-dimensional Navier-Stokes solver ; has been applied to a number of prob- 
lems of interest to GE Aircraft Engines Co, This validation exercise had two purposes The first goal 
was to evaluate the accuracy, ease of use and cost of computing the complex internal, external ant cou- 
pled flow problems with the FL3D solver. The second motive was to apply the solver in the preliminary 
mode and evaluate its adaptability to propulsion component design. A limited number of examples have 
been solved and the solutions are compared with available experimental data. The advantages and limi- 
tations with the FL3D solver are discussed Experience gained with the solver and its use in the design 
environment are welt led. 


INTRODUCTION 

The propulsion system design involves component flows 
ihai are internal, external and regions where these two 
t>pes of flows interact. In addition solutions are sought 
*hcre complex configurations such 2 l> center- bodies, mixer 
chutes, primary and secondary flow interaction, plume 
interference, boundary layer and shearlayer developments 
eic are common elements. Often in the design process, a 
number of geometric and fluid dynamic parameters art var- 
ied to understand the design performance under design and 
oil design conditions. With intuitive knowledge and under- 
standing of the flow the designers may choose to introduce 
novel changes and this may require computations to be per- 
formed with rapid turn around yielding plausible directions 
to Lhe design efforts, for example, film cooling with blow- 
ing on an airfoil at select locations. For a CFD solver lo be 
used in a general design environment, one would like it ic 
be robust, accurate, flexible and be able to handle not only 
complex configurations but also changes to the original ge- 
ometry the design process may suggest in a umely fashion 
FL3D is a robust compressible flow Navier-Stokes 
solver developed at NASA Ames Research Center The 
solver has been validated for external high-speed flows GE 
Aircraft Engines recently acquired this solver. The solver 
has been undergoing validation and call bran on tests on a 
number of internal, external and mixed flows of interest. 
Some of these problems have either experimental or other 
numerical solutions. These problems exhibit the real world 
complexiues but are easier to model. Other problems are of 
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design and explorauon type and were chosen to explore the 
user fnendJy and stahiiity aspects of the FL3D solver. 

The primary objective of this paper is to catalogue 
the experiences and understanding gained with the FL3D 
solver. As such the qualitative aspects of the solution char- 
acteristics for each case are discussed. The successes as 
well as the shortcomings of the solver are outlined. 

CFD SOLVER 

All the CFD solutions to be presented in this paper were 
obtained with the FL3D ~ 3 solver. This is an implicit, 
full Navier-Stokes solver and uses Roe upwind, TVD. 
flux-difference spbtung technique to represent the inviscid 
fluxes arid a central-difference approximauon to represent 
the viscous fluxes It is second-order accurate in space and 
first order accurate in ume. With the use of local time step 
to converge, faster to steady state, the ume ^curacy is not 
maintained. Though ume accuracy is not maintained, the 
scheme is a time marching hyperbolic method wherein the 
steady state solution is obtained oy marching from a pre- 
scribed solution The solver is capable of handling single, 
zonal or muluple grids with and without hole regions. A 
single block version with blank region (Iblanks) topology 
was used for ail the present computauors. A large num- 
ber of boundary condiuons are available and are specified 
through input by the user. These boundary conditions can 
accommodate many of the real world problems. The FL3D 
code is configured to require minimal memory and has been 
ported to Workstation platforms, such as HP workstation. 

The convergence to steady state is problem dependent 
and for simple 3-D problems like nozzles, steady state so- 
lutions have been obtained within 1000 iterations. For com- 
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plex problems, more than 6000 iterations in time were re- 
quired to reach a steady state solution. The limiting mech- 
anism seems to be the evolution of separated flow domain 
and /or the shock structures in the computauonal domain. 
The convergence rates of the solver seems to be slow while 
solving the incompressible flows (Mach number less than 

0.3 in the freestream) while maintaining good accuracy of 
the solution. 

GRID GENERATION 

A number of grid generation tools were used to gen- 
erate the grids. These tools varied from simple algebraic 
grid generators to commercially available software pack- 
ages like IDEAS and GRIDGEN3D. For complex problems 
grid generation has been the most time consuming of tasks. 
In addition to the initial grid generation, adaptive regrid- 
ding has been utilized in a limited number of problems to 
improve solution accuracy . SAGE 4 , an algebraic solution 
adaptive grid redistribution scheme, was used in the present 
work to generate adaptive grids. 

TURBULENCE MODEL 

The FL3D solver contained the original Baldwm-Lomax 
algebraic turbulence model. This model was augmented 
with a free shear-layer model This is a simple mixing 
length model and it is similar to Prandtl’s mixing length 
model for free shear layer flows. The model was developed 
with the aid of extensive data base on supersonic plumes 
measured at Florida State university 5 6 7 . In particular 

= Ci(p) * (L 2 ) * U ) (1) 

where ; 

Q = 5.9725, 

p = density, 

= vorticity, 

L = 0.60 * ( 1 + , pL) . 

V Pm f 

and 

j ] X X T 

Xa, — C 2 * r^, * tan '(— * (- — (-)cutoff + 

C 2 = 6.3662 

It should be pointed out tha, (f) CIlfo// is problem de- 
pendent . For the test cases that were studied in the present 
work (f)cutoj/ was assumed to be 20.0 . It appears to us 
that this needs to be investigated further to determine the 
validity of this model. 

The function X tq is defined in an analogous fashion sim- 
ilar to mixing lengths in boundary layer estimates. The con- 
stants and functional relationships are developed from the 
data base from Florida State University. This model was 


used to simulate the supersonic plume characteristics and 
yielded satisfactory agreement with measurements. The 
purpose of this simple model is to gain some understand- 
ing of mixing in high speed high temperature (1000 degrees 
F and up) plumes and determine dominant factors influ- 
encing the mixing. The desire is to provide data base for 
higher order turbulence models that can better describe the 
characteristics of supersonic jets. These algebraic models 
available in FL3D solver are considered ad hoc models at 
best. The limited number of turbulent solutions obtained 
with this model have proven to be cost effective and suffi- 
ciently accurate for design purposes. 

RESULTS AND DISCUSSION 

A number of internal and external flows and plume- 
interaction flow problems were computed. Summary of the 
results are given below. Computed results are compared 
with experimental data whenever available. 

1. 3D • 29° Vectored Nozzle : 

The objective is to validate the FL3D solver. In this 3-D 
internal flow problem the cross section of the nozzle tran- 
sitions from circular cross section to rectangular section at 
the nozzle exit. Also, the nozzle design is such that the out- 
flow is turned 29° to the incoming flow to vector the flow 
in a desirable direction. The flow is fully 3-D. The nozzle 
geometry , grid and mach number contours are presented in 
Figure 1 . Comparison of wall pressure distribution between 
the computed solution and the experimental measurement, 
shown in Figure 2, is very good. The mach number con- 
tours along the symmetry plane shows the flow to be very 
three dimensional. The predictions were able to capture the 
shocks ano their reflections in the aft section of the nozzle. 

Detailed look at the solution near tie throat section of the 
nozzle revealed the shock-boundary layer interaction and 
flow separation. In the design of nozzles, it is very useful 
to know the existence and extent of separated flow regions. 
The computed solution was able to predict a small recircula- 
tion zone at he throat of the top wall as expected. The good 
agreement ot the measured data and the computations are 
presented in Figure 2. Several significant features of the so- 
lution predicted by FL3D solver arc: 1) its ability tocapture 
the shocks and their subsequent reflections off of the walls, 
2) its ability to simulate shock boundary layer inter-action 
and its ability to predict the separation and 3) reauachment 
of the flow under the influence of changing pressure gradi- 
ents m the flow domain. This example is presented here as 
a validation of the code’s predictive capabilities. 

2. 3D Multi Stream S-Duct Flow • 

The objective in computing this problem was to evaluate 
the code’s capability to simulate, a) mixing of hot and cold 
streams, b) compression comer effects, and c) expansion 
comer and curvature effects. The geometry for this problem 
was selected carefully to include all the features to simulate 
the above mentioned complexities. At the inlet an axisym- 
mctric cold-stream at 880° R and a hot core flow at 1995° 

R were introduced. These two flows are separated at inlet 
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by a wall. The schematic of the flow features, the geometry 
and a representation of the grid used are shown in Figure 3. 
The gnd is dense at the region where the cold and hot flows 
meet in the inlet region. Both the flows at entrance are sub- 
sonic. The mixed flow accelerates to Mach number 1.0 as 
it approaches the throat. Bevond the throat region the flow 
supersonic flow as it passes through the noz- 
zle. The cross sectional gnd, total temperature distribution, 
Mach contours along the axial plane and at selected cross 
sections are presented in Figure 4, The Mach number distri- 
bution shows a large acceleration on the expansion comer 
and deceleration at the bottom comer. The total tempera- 
ture distnbution shows the develoomem of the shear Laver 
and the mixing of the hot and cold flows in this region. T»e 
sharp curvature changes in the geometry gave rise to the 
formation of two vortices. These vortices in turn promoted 
further mixing of the two streams in the ducL The total 
temperature contours presented in Figure 4, clearly shows 
the dramatic effect of these vortices. Beyond the throat a 
divergent nozzle was attached to the geometry to further 
examine the capabilities of the code, FL3D, to simulate the 
supersonic flow. Details of the calculations showed the flow 
separation and reatiachment at the top wall of the duct The 
recirculation regions created by the compression comer at 
the bottom wall are also captured satisfactorily by the pre- 
dictions. 

In summary the qualitative results presented here demon- 
strate the ability of the FL3D code to predict the mixing of 
hot and cold flows coupled with the influence of the sharp 
curvature changes in die geometry. It shows further evi- 
dence of the code’s ability to satisfactorily simulate the in- 
fluence of expansion corners in supersonic flow regimes. 

3. 2-D Splitter Plate flow: 

This problem was selected to demonstrate the influence 
of geometry changes on the solution. The georr’try con- 
sists of channel formed by a 40 inch long top and bouom 
walls. A splirer plate of 8 inches long is introduced to split 
this channel in two separate channels at the entrance. In 
the top channel a cold subsonic stream was introduced. In 
the bottom channel a hot choked flow was introduced. The 
static temperature ratio of the two flows (Hot/Cold) is 1.1 . 

In this exercise two problems are considered. For both 
the problems the bottom wall and the splitter plate geome- 
tries are held at fixed positions. In the first problem the 
top wall is held parallel to the bottom wall. Whereas in the 
second problem the top wall was diverging loan overall in- 
crease of 5 per cent in the exit area starting at 20.0 inches 
from the inlet. The hot flow between the splitter and the 
bouom wall is expanded (transonic flow mach 1.2 at the 
end of splitter plate ) and is interacting with the coid flow 
from the top channel. The development of the boundary 
layer on the top wall and the mixing layer from the end of 
the splitter plate are captured properly by the code. 

The splitter plate thickness and small changes in the an- 
gle of the top wall had significant effect on the flow devel- 
opment. The computational gnd and the temperature con- 


tours from the computed solutions arc shown in the top of 
Figure 5. The bottom two figures compare the estimated 
temperature distribution on the top wall. The small changes 
in geometry had significant effect on the solution and this is 
demonstrated by changing the top channel (secondary flow) 
cross-section slightly. This in turn caused large changes in 
the predicted pressure. The primary(bottom clumnel)/core 
flow can be allowed to reach either subsonic or supersonic 
conditions by the small changes introduced in the secondary 
(top channel) flow region. The differences in the computa- 
tions and the measurements are due to the uncertainty in the 
geometric data provided by the experimental groups. Ge- 
ometry for the first 20 inches are known accurately and the 
rest of the geometry was not The disagreement of die com- 
putations with data for this example may be parJy due to the 
turbulence model and its inability to simulate supersonic 
mixing satisfactorily. 

4. Supersonic nozzle flow with center body: 

The aim here is to introduce a contoured center body in a 
supersonic nozzle and evaluate the solution. The geometry 
considered is an axisymmetric nozzle with center body and 
the flow is choked at the throat and the flow is supersonic 
at the exit. Figure 6 shows the grid and the Mach number 
contours from the predicted solution. The choking of the 
flow at the throat and the compression produced along the 
compression region of the center body due to the change in 
curvature is clearly seen. A shock wave has developed at 
this region on the center-body (where the curvature of the 
geometry is changing) and is subsequently reflected off of 
the nozzle wall. The effects of the curvature on supersonic 
nozzle flows were caouired well by the code with no sur- 
prises. 

5. The effect of External Pressure: 

The intent of this problem is to investigate the effect 
ot the over-expansion of the supersonic nozzle flow under 
simulated freestream flow conditions. An external (out- 
side wall) wail was added to the nozzle geometry shown 
in the previous problem. A high subsonic free stream flow 
of Mach number 0.9 is included in the problem definition. 
With this external flow condition the nozzle exit flow was 
over expanded. The combined grid and the FL3D estimated 
mach number contours are shown in Figure 7. The interac- 
tion of these flows (nozzle flow and external flow) created 
an extremely complex flowfield as presented in the mach 
number distnbution shown here. The external wall of the 
nozzle is contoured to accelerate the freestream flow on it 
as a result of this flow acceleration on the nozzle external 
wall Mach number exceeds 1 .0 on it, thus reducing the pres- 
sure on the wall. The subsequent recompression of the flow 
to regain the static pressure balance gives rise to the shock; 
a lambda shock; on the nozzle wall. The shock-boundary 
layer interaction, coupled with the over expanded flow led 
to tfie complicated shock pattern inside the nozzle as shown 
in the figure. These interactions separated the flow from the 
nozzle wall both inside as veil as out side. A plot of the ve- 
locity vectors, Figure 8, at the trailing edge of the nozzle 
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clearly indicate the flow reversals, recirculations, and sep- 
aration due to the shock boundary layer interactions in an 
over expanded supersonic flow under simulated high sub- 
sonic freestream flow. 

6. Axisymmetric Turbulent Plumes: 

Plumes are an integral part of propulsion system. They 
can occur either external or internal as in the case of sin- 
gle sided expansion ramp, one sided evolution of plumes 
determine the thrust efficiency. To validate the CFD solver 
for plumes axisymmetric under-expanded and ideally ex- 
panded plumes are considered. Detailed measurements of 
the flow characteristics from such supersonic jets are made 
at Florida State University (1991-1992 Mechanical Engi- 
neering Department/ GE - FSU Cooperative program). A 
simple mixing length turbulence model was devised to sim- 
ulate the turbulent jets. Calculaaons are made with this 
model to predict the characteristics of the under-expanded 
jets and fully expanded jets. The predicted temperature 
contours for these jets are shown in Figure 9. The shock 
patterns shown for the under-expanded jets are indicative 
of the plausibility of the solution. The shock free nature 
of the fully expanded jets is also presented in this figure. 
A comparison of Mach number distribution on the center- 
line of the jet with measured data is presented in Figure 10. 
The agreement between the measurements and predictions 
for the fully expanded (shock free conditions) jets shown 
in figure is very good. For underexpanded jets, the noz- 
zle pressure is higher than the ambient pressure. As such 
the jet expands to reach a pressure balance with that of the 
ambient pressure; thus giving rise to higher Mach number 
down stream of the nozzle exit In this process the jet over 
expands and tnes to achieve pressure balance via compres- 
sion (oblique shock formation) wave. Tnis procession of 
expansion and compression, after the nozzle exit contin- 
ues till the pressure in the jet reaches that of the ambient 
pressure. The ability of FL3D to predict such flow proper- 
ties with proper frequency ( the number of expansions and 
compression cells ) and amplitude (magnitude changes) of 
the mach number is shown in the figure. Comparison of the 
predictions with experimental data presented in the figure 
shows satisfactory agreement. 

These comparisons demonstrate the ability of the FL3D 
to accurately predict a variety of complex flow fields from 
propulsion systems 

Conclusions 

The predictions and their comparisons with experimen- 
tal data; whenever available; demonstrate that the code 


FL3D can capture the salient characteristics of super- 
sonic/transonic nozzle flows. Curvature effects of the noz- 
zle geometry on the flow fields are predicted accurately. Its 
ability to estimate the properties in a shearlayer is demon- 
strated in the splitter-plate problem. The computational ex- 
periments performed further show the capabilities of the 
code to predict separation, recirculation and flow reversals 
in a complex flow field. Its ability to capture shocks has 
been demonstrated amply through the exercises presented 
here. The application of code to predict the supersonic jet 
flow fields has been demonstrated in this work. Our ex- 
perience shows that FL3D can simulate the flow fields for 
propulsion system? adequately. 
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Figure 1 . 29° vectored nozzle geometry, grid and 
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Figure 2. Comparison of wall pressure distribution 





Problem: Investigate mixing of subsonic primary (hot) flow with 
subsonic secondary (cold) flow exiting as supersonic 

Nature of the Flow: Transonic, mixing of hot core and cold fan 
flow with sharp curvature 
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Figure 3. S-Duct grid and flow features description 
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Figure 4. S-Duct Flow characteristics 





Problem: Investigate the mixing of supersonic flow with a subsonic 

film at one wall 
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Figure 5. Mixing characteristics of supersonic 
and subsonic flows 




Problem: Supersonic flow solution with center body 

Nature of the Flow: Transonic 
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Figure 6. Supersonic nozzle flow and geometry 



Problem: Investigate the flow properties from an over-expanded 

supersonic nozzle 

Nature of the Flow: Over expanded supersonic flow at 0.9 external 
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Figure 7. Over expanded nozzle flow 
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Figure 9: Supersonic plume characteristics 




Problem: Estimate the plume characteristics of a supersonic jet 

Nature of the Flow: Fully expanded and under expanded plumes 
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ABSTRACT 


A two-dimensional, convergent-divergent (2DCD) exhaust system has been analyzed with the FL3D 
Navier-Stokes how solver. The nozzle incorporated a chute suppressor system at the exit of the core 
flowpath for noise abatement at takeoff. Results show an extremely three dimensional dowheld with 
very complex how structures. The computed solution was compared to limited experimental test data. In 
general ' the predicted pressure distributions matched well with test data. The predicted flow field captured 
the important flow- features observed during testing and from shadowgraph pictures. The computational 
results have enhanced the understanding of this type of complex nozzle flow and can be used in improving 
its design. 


INTRODUCTION 

Recent studies have predicted a 400 percent increase in 
commercial transpacific air travel and a doubling of transat- 
lantic tops by the year 2000. These studies have shown 
a sizable market potential for a civil aircraft which would 
significantly reduce the ume to fly these transcontinental 
routes. An aircraft that cruises at over twice the speed of 
sound will cut the ume to fly from Los Angeles to Tokyo 
irom 10 hours to about 4 hours. In a world that is truly 
becoming "globalized", this ume savings is very attractive. 

One of the challenges to high speed aircraft systems is 
designing an engine system which meets or exceeds civil 
noise requirements. During takeoff and climb, the domi- 
nant noise source emanates from the engines. Typical mil- 
itary engines for supersonic flight incorporate low bypass, 
high specific thrust turbofan systems. As such, these sys- 
tems generate high exit velociues to meet thrust require- 
ment. These high exit velocities in turn generate high lev - 
els or noise which are unacceptable in a civil aircraft. In 
order to lower the noise levels resulting from these engine 
systems, exhaust nozzles which entrain ambient air into the 
nozzle are a prime candidate for this aircraft 1 ” 10 . Ihc en- 
trainment and mixing of this ambient air with the engine 
core flow lowers the exit velocity while maintaining the 
thrust at required levels. 

This paper presents a computational study performed to 
evaluate one of these mixer nozzle candidates. The FL3D 
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flow solver from NASA Ames Research Center was used 
for this analysis. As will be described in more detail later 
in the paper, FL3D is a three-dimensional (3D) full Navier- 
Stokes computational fluid dynamics (CFD) solver. The 
main purpose of this study was to better understand the 
flow physics of this nozzle design while validating FL3D 
for these type of systems. The results of the analysis are 
compared to a limited amount of test data and several inter- 
esung charactensucs of the nozzle flowfieid are discussed. 

2DCD NOZZLE GEOMETRY 

The nozzle analyzed in this study is shown schematically 
in Figure 1. It is a two-dimensional (2D) wedge nozzle 
with a lobed mixer located in the primary stream. A similar 
nozzle was designed and tested by Harrington, et al. 1 The 
unique feature of this nozzle is the addition of a secondary 
or bypass stream to further attenuate jet noise. 

The mixer rack consists of five mixer chutes located at 
equal spanwise distances apart As can be seen in Figure 1 , 
the chute entrance angle is relatively steep (55 degrees). 
This was dictated by the need to stow these chutes inside 
the wedge at cruise conditions. The orimary flow flowpath 
through the mixer is convergent, resulting in the primary 
stream throat being located at the mixer exit. At the exit of 
the mixer (i.c. the mixing plane) the secondary duct area is 
1.5 times the primary duct area. 

The nozzle has a straight wedge with a 15 degree half 
angle. The sidewalls for the secondary duct terminate at the 
secondary nozzle exit. The sidewalls for the primary duct 
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terminate ai the mixing plane. The secondary flowpath is a 
convergent nozzle with the throat at the exit plane. 

The spanwise computational domain for this analysis 
(shown in Figure 2) consisted of a section of the total nozzle 
geometry extending from the centerline of a primary chute 
to the centerline of a secondary chute The top- to- bottom 
computational domain extended from a farfield freestream 
boundary to the nozzle centerline. 


GEOMETRY MODELING USING I DEAS 

Since the nozzle geometry is highly three-dimensional 
an advanced CAE package called I-DEAS (Integrated De- 
sign Engineering and Analysis Software) was used to create 
the geometry model and the computational grid for the flow 
analysis. The nozzle geometry is generated using a build- 
ing block approach provided in the I-DEAS solid modeling 
environment. At first the chute geometry is built using the 
cross sectional definition and invoking the ’sinking opera- 
tion ’ . Then a thickness is added to the inner chute surface to 
generate the actual chute geometry with variable thickness 
from top to the bottom of the chute. Then the single chute 
is replicated to give a series of chutes that represent the true 
chute configuration in the nozzle. 

The nozzle flow domain without the chutes being two 
dimensional, was created using the outline of the inner and 
outer boundaries only. The chute is then located in this 
domain in its proper orientauon and location to generate 
the divided flow path for the core and the fan flow. The 
most important information in this operation is the three di- 
mensional intersection of the plug boundary with the chute 
/ boundary at its bottom. This information is automatically 
y\ obtained from I-DEAS when the chute is located or. to the 
' plug surface. The 3-D surface information and the outer 
flow boundaries are then transferred directly to the grid gen 
erator for creating suitable gnd for the flow analysis. 

COMPUTATIONAL GRID 

Since the geometry is complex with immersed bodies 
and wall boundaries in the flow field, a suitable block struc- 
ture was designed to accommodate these complexiues in 
the flow domain. The block structure is also designed such 
that a better control of grid in terms of quality as well as 
ts . quantity could be achieved in the flow domain. A schematic 
V of the grid blocks generated is presented in Figure 3 . The 
blocks are generated using the block structure opuon avail- 
able in I-DEAS gnd generation module. This block struc- 
ture algorithm through data base management system keeps 
track of the block interface connectivity between two differ- 
ent blocks as well as the surface interface connectivity of a 
single block. This feature automatically allows the same 
number of grid points and the function continuity along 
the interface between the blocks. The grid is generated m 
these individual blocks first and then merged to form a sin- 
gle structured grid. Figure 4 shows the computational grid 
in the chute-symmetry plane. A finer grid distribution is 


provided where ever there is a solid boundary to resolve 
the viscous layer near the wall. The grid was clustered in 
the stream wise direction where ever large gradients are ex- 
pected in the flow field such as the region close to the exit 
plane of the primary nozzle. 

RESULTS AND DISSCUSSION 

An algebraic turbulence model based on a Baldwin- 
Lomax formulation has been used with the FL3D flow 
solver for the results presented in this paper. The length 
scales are deduced from the grid normal distances from the 
wall, and the velocity scales are estimated from the venic- 
ity distribution in the flowfield. The results presented were 
computed both with laminar as well as with turbulent vis- 
cosity. The influence of turbulence on the pressure distribu- 
tion on various surfaces of the geometry is explored. The 
changes in the flowfield from the core center plane to the 
chute center plane are presented. The interaction of the core 
flow with the secondary flow and its impact on the separa- 
tion characteristics of the flow on the centerbody is inves- 
tigated. The estimated pressure distributions are compared 
with measured data where available. Qualitative compar- 
isons of computed flow characteristics with shadowgraph 
measurements are also shown in this paper. 

FLOW ORIENTATION 

The flow orientation presented earlier in the paper is 
briefly re-stated here. The computational domain consists 
of three different flows. They are : (1) Primary flow, (2) 
Secondary Flow and (3) Freestream flow. These flows are 
separated by solid walls up to the nozzle exits. Downstream 
of the nozzle exits, they are free to interact with each other. 
The physical orientation of the 2DCD nozzle geometry and 
the flow are shown in Figure 1. The flow conditions used 
i in this study are: 

~\ Primary Pressure Ratio (PT p „„,/P«mt>) = 3.2 
Secondary Pressure Ratio (PT«ec/P«mb) = 2.2 
Freestream Mach number = 0 2 

The primary flow accelerates to supersonic speeds as 
it negotiates through the passage between the chutes. As 
such, the pressure of the primary flow at the mixer chute 
exit is dictated by the area ratio of the primary chute flow- 
path. However, at the end of the mixer chute, the flow exits 
to ambient conditions and must adjust accordingly. 

The secondary flow exits the nozzle at supersonic speeds 
just upstream of the mixer chute leading edge. The flow is 
then split into two regimes. Some of the flow is free to flow 
down the chutes and mix with the primary flow in the adja- 
cent mixer chutes. The remaining flow continues over the 
top of the primary chutes where it forms a stratified layer 
over the primary flow. The freestream flow is able to in- 
teract with the secondary and primary flows downstream of 
the secondary nozzle exiL 




AIAA Paper No. 94-0019 


lid to a large scale flow reparation on the wedge surface, 
and separation of the secondary flow in the mixer chutes. 
The turbulent flow results agreed well with available sur- 
face static pressure measurements and with shadowgraph 
pictures of the flowfield. The CFD results have provided a 
good understanding of the flow phenomenon in these type 
of nozzles and will aid in improving its design. 
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Figure 1- Schematic Ol 2DCD Nozzle Geometry 


Secondary Chute 



Pnmary 

Chute 




} 


Secondary 

Chute 


Figure 2. Schematic of the computational domain at the mixer exit. 



Figure 3. Block structure representation for the computational grid. 
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Figure 8. Computed velocity vector direction - core-symmtery plane. 
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Figure 9. Computed velocity vector direction - chute-symmtery plane 
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Figure 11. Particle traces from fan flow region alone. 
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Figure 12. Pressure distribution comparison along the wedge (ramp) wall 
- laminar computations. 
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Figure 13. Pressure distribution comparison along wedge wall 
- turbulent computations. 
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Figure 14. Pressure distribution comparison along the wedge (ramp) wall 
- laminar computations. 
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Figure 15. Pressure distribution comparison along wedge wall 
- turbulent computations. 
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Chute Trailing Edge Pressure Distribution 
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Figure 16. Pressure distribution comparison along the chute trailing edge 
- turbulent computations. 
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Figure 17. Pressure distribution comparison along the chute center-line 
- turbulent computations. 
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Figure 1 8. Comparison of the computed Mach contours on the core-symmetry 
plane (left) and the experimental shadowgraph (right). 
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Figure 19. A composite view from the computations : surface pressure, 
particle traces and Mach contours on axial plane downstream of the chute 
exit. 
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Abstract. An initial experimental/numericai investigation has been conducted to gain a better 
understanding of the multi-dimensional flow phenomena inside pulse facilities and the influence 
of these phenomena on test conditions and test times. Experimental data from the NASA Ames 
electric-arc driven shock tube facility (from cold driver shots) is compared to time-dependent 
axisymmetric numerical simulations of the complete facility. These comparisons help establish 
the numerical modelling requirements for simulating shock tube flow and help validate the com- 
putations. The numerical simulations are used to study the interaction between the reflected 
shock wave and the side wall boundary layer and the resulting shock bifurcation. Of particular 
interest is the effect of the bifurcated shock structure on the driver/driven gas interface. The 
computations demonstrate how this shock structure introduces a mechanism for the driver ga: 
to contaminate the the stagnation region thereby reducing the duration of the test time. The 
simulations incorporate finite-rate chemistry, a moving mesh and laminar viscosity. 

Key words: Shock Tube, Computation, Shock/Boundary Layer Interaction 

1. Introduction 

Knowing the length of the test time and the state of the test gas provided by shock tubes or 
shock tunnels is critical to interpreting data obtained from these facilities. Unfortunately, it is 
well established that the test time achieved in these facilities is usually significantly less than 
ideal theory predicts. Many investigations have been carried out to understand and quantify the 
physical mechanisms which cause shortened test times. These mechanisms include deformation of 
the contact discontinuity caused by the diaphragm rupture process, mass transfer of the driven gas 
into the boundary layer, contact discontinuity instabilities, and shock/boundary layer interaction 
after reflection of the incident shock off the end wall. 

It is not clear which of the mechanisms mentioned above limits test times most and it may be 
that the dominant mechanism varies with the experimental facility or the run conditions. There 
is a large amount of evidence, however, that suggests that the reflected-shock/boundary layer 
interaction is often a major contributor to the contamination process. It is well established that 
under many conditions the reflected shock will interact with the boundary layer causing it to 
bifurcate near the wall (see Fig. 1). Mark 1958) developed a model to predict the characteristics 
of the bifurcation and the conditions under which it will occur. He showed that the flow in 
the energy deficient boundary layer has a stagnation pressure that is less than the stagnation 
pressure behind the normal reflected shock and is prevented from passing under the reflected 
shock. Instead, it separates and collects in a bubble of gas next to the wall. A conse uence of 
the shock bifurcation is a jet of gas near the shock tube wall as depicted in Fig. 1. The jet is 
present because the gas which passes through the oblique shocks a* the foot of the bifurcated 
shock retains a higher velocity than the gas which encounters the normal shock. Particularly clear 
experimental evidence o* this phenomenon can be seen in the recent color schlieren photographs 
of Kleine et al. (1991) which include features such as the rolling up of the wall jet as it encounters 
the end wall. Davies (1966, 1967, 1969) used Mark’s model to show that the wall jet provides a 
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mechanism for contamination of the stagnation region by propelling cold driver gas toward the 
end wall and into the driven gas. He and others such as Bull and Edwards (1968) have done 
experiments which measured the time of arrival of the cold driver gas through this mechanism. 




Fig. 1 . Schematic diagram of reflected shock/boundary Fig. 2. Schematic diagram of pressure trace under 

layer interaction bifurcated shock 


The reflected-shock/boundary layer interaction described above can explain much of what is 
observed experimentally in shock tubes. However, more complicated flow structures such as a 
pseudo-shock (or shock train) can develop when the shcck/boundary layer interaction is strong. 
Matsuo et al (1975), Strehlow and Cohen (1959), and Brossard and Charpentier (1985) all show 
schlieren photographs of the formation of multiple shocks after the reflection of the incident shock. 
The effect of these multiple shocks on driver gas contamination has not been studied. 

In the absence of optical data, the presence of shock bifurcation can be inferred from side- wall 
pressure measurements. Figure 2 shows that the passage of the reflected shock is marked by a 
two-step pressure rise follow by a pressure overshoot as described by Sanderson (1969). It can 
also be deduced by noting the change of the reflected shock speed compared to inviscid theory 
(Mark 1958; Strehlow and Cohen 1959; Matsuo et al. 1975). 

There have been several recent examples of computations of reflected-shock/boundary layer 
interaction at the end of a shock tube with the computational domain limited to the stagnation 
region (e.g. Kleine et al. 1991 and Yakano 1991). The present authors are unaware, however, of 
any computations which have looked at the contamination of the stagnation region with driver 
gas through the wall jet mechanism proposed by Davies. One way to accomplish this is to begin a 
simulation of a shock tube at the diaphragm rupture allowing the position of the contact discon- 
tinuity and the boundary layer development to be computed. This, in turn, makes it possible for 
phenomena such as the reflected shock/boundary layer interaction and the reflected-shock/contact 
discontinuity interaction to be investigated numerically. This approach is adopted herein by com- 
puting the time-dependent flow inside the NASA Ames electric-arc driven shock tube (cold driver 
shots without the arc-driver). The simulations assume that the contact discontinuity is planar at 
diaphragm rupture and that boundary layer is laminar. Experimental data in the form of wall 
static pressure traces and heat transfer was gathered to help guide the numerical modelling and 
validate the simulations. The present work is a continuation of research reported in Wilson et al. 
(1993) and more details can be found there. 
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2. Experimental Facility 

The NASA Ames electric-arc driven shock tube facility has several possible configurations and has 
a large hypervelocity operating range using its arc-driver (Sharma and Park 1990). However, this 
work only considers experiments using a cold helium driver with nitrogen in the driven section. 
These shots were made with a cylindrical driver of .86 m (2.8 ft) in length and 10 cm (3.93 inches) 
in diameter. The driven section was 4.22 m (13.85 ft) long with the same diameter as the driver. 

A single self-break diaphragm separated the driver and driven gases. 

For the current experiments the instrumentation consisted of pr^ure transducers flush 
mounted at three fixed positions on the shock tube walls, at 2, 6, and 24 inch > from the end 
wall. There were also two pressure transducers flush mounted on the end wall, one at the center 
and one 1/4 inch from the tube side wall. In addition, the end wall was modified so that it could 
be moved forward and backward at one inch increments. This made it possible to collect data at 
variable distances from the end wall (i.e. 1, 2, and 3 inches from the end wail for the first side 
wall gauge and, thus, 5, 6, and 7 inches from the end wall for the second gauge). The pressure 
transducers were PCB Piezotronics, Inc. Model 113A21 with a circular surface area .218 inches 
in diameter and a rise time of 1 microsecond. Heat- transfer data has also been collected at the 
same side wall locations at instrument locations opposite of the pressure transducers. 

3. Numerical Method and Gas Model 

The gas dynamic equations for the axisymmetric simulations are solved by using an explicit finite- 
volume form of the Harten-Yee upwind TVD scheme (Yee 1989). The ga model includes the three 
major species present in the shock tube for the present experiments (JV^, N , and He) and accounts 
for finite-rate chemical processes. A separate equation for vibrational energy is included so that 
vibrational nonequilibrium effects can be assessed- The present work, however, enforces thermal 
equilibrium. The numerical method is essentially an extension of the quasi-one-dimensional work 
in Wilson (1992). The full Navier-Stokes viscous terms are included. 

Mesh points are clustered at the contact discontinuity to minimize numerical diffusion and 
are convected with the gas interface as it travels down the driven tube. This approach has the 
additional benefit of compressing all of the cells associated with the driven tube into the end wall 
region of the shock tube as the driven gas is compressed thereby providing a fine axial mesh during 
the shock reflection. The grid is also clustered around the incident shock so that the resolution 
is relatively high where the boundary layer is initially formed behind the shock. Points are also 
concentrated near the wall to resolve the boundary layer. Because the Euler terms are treated 
explicitly, the computations are advanced at a CFL number less than 1 based on the inviscid 
gas dynamics. To avoid the more limiting time step dictated by the viscous terms, the thin-layer 
viscous terms are treated implicitly (note that all the viscous terms are included explicitly). It 
is believed that the time accuracy of the solution is not significantly effected by this approach 
and without it the simulation becomes impractical. The source terms representing the finite-rate 
chemical kinetics and vibrational relaxation art also treated implicitly. This implicit formulation 
reduces the formal temporal accuracy to first order. 

5. Results and Discussion 

The results presented in this short paper represent a single test condition with a nominal driver 
pressure of 4.83 MPa (700 psi) and a driven pressure of 20 torr nitrogen with both sections at 
ambient temperature (295 K). Three shots are reported with the end wall position varied for 
each shot so that there is data 1, 2, and 3 inches from the end wall (these shall be referred to 
as the 1 inch position, the 2 inch position, etc.). Because of the self-break single diaphragm, a 
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precise repeatability of shock speed for the three shots was not possible. The shock speed for the 
1, 2, and 3 inch positions were 1984 m/sec, 2005 m/sec. and 1953 m/sec, respectively. Figure 3 
shows the pressure trace for the 2 inch position. Several of the major events are labeled. These 
include the passage of the incident and reflected shocks, the arrival of the rarefaction, and the 
presence of waves reflected off the contact discontinuity. The time axis is adjusted so that zero 
time corresponds to the shock arrival at the end wall. This allows a meaningful comparison of 


traces with the end wall in different positions. 



Fig. 3. Experimental and computed pressure traces 2 
inches from the driven tube end wall 



Fig. 4. Experimental and computed pressure traces 
1, 2, and 3 inches from the driven tube end wall 


The axisymmetric simulation used an initial driver pressure 1 1 % higher than the experimen- 
tally reported value in order for the computed incident shock speed to match the experimental 
one. The reason for this discrepancy is not known. All other initial conditions matched the experi- 
mentally reported ones. A cold wall boundary condition of 295 K was enforced. The computational 
mesh contained 800 points along the length of the tube (400 each in the driver and driven sections) 
and 112 points between the tube centerline and an outer wall. The points were exponentially clus- 
tered near the wall with a minimum spacing of .150 mm at the wall for the first few meters of the 
driven tube and was ramped down to a constant .015 mm over the last .75 meters of the tube. 
This approach eased the computational cost by allowing a 1 rger time step early in the solution. 
As mentioned before, the solutions assumed an initially planar contact discontinuity and laminar 
viscosity. Verification that the reflected shock/boundary layer interaction has all of the features 
depicted in Fig. 1 is presented by Wilson et al. (1993). 

A composite of experimental and computed pressure traces in the end wall region for early 
times after the shock reflection is presented in Fig. 4. It is seen that the general features of the 
experimental and computed pressure traces are quite similar. Evidence of the shock/boundary 
layer interaction in the pressure traces, as depicted in Fig. 2, is clearly seen. The computation 
and experiment are in good agreement for the trace 1 inch from the end wall, fair agreement 
at 2 inch position, and by the 3 inch position it is clear that discrepancies are growing. Further 
deterioration of the agreement is seen in the traces farther back from the end wall (not shown 
here). While it appears that the growth of the bifurcated shock structure is being captured fairly 
well by the simulations, the computed speed of the bifurcated shock is too fast and the predicted 
pressures under the separated flow region differ. Even with the differences, the general features of 
the flow appear accurate enough that a qualitative investigation of the interaction of the reflected 
shock with the contact discontinuity is deemed worthwhile. 
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Fig. 5. Computed temperature contours showing reflected shock/contact discontinuity interaction 


Figure 5 contains temperature contours after the reflected shock has interacted with the 
contact discontinuity. The figure shows that the bifurcated shock structure deforms the contact 
discontinuity near the wail. Hot driven gas in the separated region under the shock bifurcation 
is carried into the driver gas while the driver gas that passed through the oblique shocks of the 
bifurcation retains a higher flow velocity and penetrates into the driven gas (i.e. it is part of the 
wall jet). There is also the additional feature of a shock reflected off the contact discontinuity due 
to the overtailored conditions. There are many nonuniformities in the stagnation region because 
of the shock/boundary layer interaction. As waves such as the shock reflected off the contact 
discontinuity interact with these non-uniformities even more complex features are formed. It is 
easy to see why pressure trace data can become quite noisy. 

There are many possible reasons for the discrepancies seen between the experimental data and 
the present computations. Experimentally, the point of transition between laminar and turbulent 
flow is not known. It is believed that the flow near the end wall is laminar (heat transfer data seems 
tc support this); however, the flow must eventually become turbulent. This may happen in t he 
separated flow region under the bifurcated shock. Data interpretation is further complicated by 
arbitrary variations in the finer features. For example, the size of the pressure overshoot associated 
with the shock bifurcation can be notably different on the same shot for two transducers mounted 
at opposite sides of the tube. Numerically, the authors have found that it very difficult to resolve 
the shock/boundary layer interaction to the point wher solutions become grid independent. 
Many grid refinement studies have shown that the wail spacing used here (.015 mm) is nearly 
sufficient;, however, these studies have also shown that the grid spacing along the tube is also 
important and that the current solutions would be helped by further grid refinement. Additionally, 
the starting assumption of an sharply defined, planar contact discontinuity causes waves which 
interact with the interface to be much more sharp (and often of higher magnitude) than those 
observed experimentally. An initially deformed interface will be tried in the future. 


5. Conclusions 

Axisymmetric simulations of the NASA Ames electric- arc driven shock tube have been done which 
include the wall boundary layer and the computations have been compared with experimental 
data. These simulations have allowed the wall jet created by the reflected shock/boundary layer 
interaction to be investigated numerically for the first time. These simulations support earlier 
analytical and experimental work which indicate that this mechanism can contribute to the re- 
duction of the usable test time by allowing the driver gas to contaminate the stagnation region. 
Before any driver gas contamination, the wall jet creates nonuniformities in th<^ stagnation region 
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which create complex flow patterns, especially when waves due to non- tailored conditions interact 
with these non-uniformities. 
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Extended Abstract 

We have used laser induced fluorescence (LIP) to characterize the free str eam of the 20 
Megawatt Aerodynamic Heating Facility at NASA/' Ames Research Center. Knowledge of the 
free stream, properties is required for a detailed understanding of the interaction between the arc 
flow and test models. Until recently the properties of the free stream have been inferred from 
measurements of bulk properties such as heat loss to cooling water in the discharge and heat 
transter to test models, combined with emission spectroscopy and computer modeling of the arc 
jet flow. He. : loss measurements in the discharge are too difficult to be carried out routinely, 
and model heat transfer measurements suffer from a lack of detailed knowledge of the chemical 
processes occurring on the surface. Because the most important species in the flow are in the 
ground electronic state, emission spectroscopy can provide no information about them. Our ex- 
perimental program is designed to measure the velocities, translational temperatures, and radial 
distributions for two of the most important species, atomic oxygen and atomic nitrogen. Com- 
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parison of the experimental results with both theory and earlier experimental results will greatly 
increase our knowledge of the free stream properties. 

The experimental approach for number density measurements is as follows. The photon 
source is a Nd: Y AG-pumped, frequency doubled dye laser producing tunable radiation in 
ultraviolet. To detect a.omic oxygen, we excite the Sp^P^p^P two-photon transition at 226 
nm and detect fluorescence from the 3p^?2 \ q 3s- transition at 845 nm. For atomic 
nitrogen detection, we plan to excite the Sp^D^p^S two-photon transition at 21 1 nm and 
detect fluorescence from the 3p^D— >3s 4 P transition at 870 nm. The fluorescence is imaged onto 
the surface of a red-sensitive photomultiplier tube. For number density measurements the in- 
tensity of the near-infrared laser-induced fluorescence (LIF) is monitored as a function of the 
wavelength of the ultraviolet laser. Using calibration procedures described previously [Phys. 
Rev. A 34, 185 (1986)1 along with absolute values of the two-photon cross section (Phys. Rev. 

A 36, 3497 (1987)] the absolute number density is determined. This measurement is repeated at 
several different radial positions (with a spatial resolution of about 1 cm) to deter min e the radial 
dependence of the number density in the flow (which has cylindrical symmetry). 

Flow velocities for atomic oxygen are determined from the Doppler shift of the 
3p^P2 i o 2p3?2 transition relative to the fixed frequency of a nitric oxide LIF transition in a 
static room temperature gas celi. We first record the unshifted position by taking several simul- 
taneous spectra of flowing O-atoms and static nitric oxide with the laser at 90 degrees to the 
flow, as shown in Figure 1. We then repeat the experiment with the laser at an angle 9 degrees 
to the flow, as show n in Figure 2. The spacing between the R j ( 1 6) line of the 7(0.0) band of 
15 n16q and the 3p^P2,l,0 2p^2 tine °f atomic oxygen is measured in both cases, and the 

difference between these spaeings, A A , is determined. The velocity', v, is then calculated using 
the following equation: 


v -= (AA/A) x (c/cos 6 ), 

where c is the speed of light. The precision of the O-atom velocity measurements is ±1 km/sec 
For atomic copper we use a similar approach, eliminating the need for a frequency standard by 
simultaneously exciting the flow with beams at 9 =90 degrees "id 0 =73 degrees. Excitation 
occurs via the ^ 3 / 2^^112 transition at 324.85 nm or the transition at 327.49 
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nm, and fluorescence is detected from either the 2 ^3/2^ 2 ' D 5/2 transition at 510.70 nm or the 

^1/2 ^ d 3/2 transition at 578.37 nm. The precision of the Cu-atom velocity measurements is 
±1.5 km/sec. 


The relative O-atom number density as a function of radial position, for expansion 

through a 61 cm diameter nozzle, is shown m Table I. 


Table I. 

Relative number 

density of 0(2p- ) as a function of radial 

. position. 

~ “ 





distance from 

relative fO] 



centerline (cm) 

(arbitrary units) 


0 


0.86 



15 


0.85 



25 


0.63 



= — = = 





In future work we will measure the relative number density at several additional radial positions 

and place the number densities on an absolute scale. 


The experimental velocities are given in Table II. 


Table EL 

Summary of experimental velocity measurements. 




nozzle 

mass-avg. 


Date 

species 

dia. 

enthalpy 

velocity 



(cm) 

(MJ/kg) 

(km/sec) 

9/16/92 

Cu 

30.5 

15 

3.0 

9/17/92 

Cu 

30.5 

15 

2.6 

9/21/92 

Cu 

30.5 

13 

3.2 

9/21/92 

Cu 

30.5 

13 

3.9 

12/17/92 

0 

30.5 

17 

2.7 

2/19/93 

0 

61.0 

15 

4.2 

2/19/93 

0 

61.0 

15 

4.0 

3/9/93 

0 

61.0 

15 

4.2 

3/10/93 

0 

61.0 

12 

5.3 

3/31/93 

0 

61.0 

10 

4.1 

3/31/93 

0 

61.0 

9 

4.4 


Ail of these velocities except the last one were measured along the centerline of the flow. The 
last velocity was measured at a point 15 cm above the centerline of the flow to test for a spatial 


dependence of the flow velocity. Within our error limits, there appears to be no spatial depen- 
dence. Because the velocities measured under (nominally) identical conditions are internally 
consistent, they can be used with confidence to test existing codes describing the arc jet flow. Ln 
the near future we will add another column to Table II containing the calculated velocities and 
briefly describe the procedure used to obtain them. 

To measure absolute number densities of 0(2p^P) we must measure the rate of col- 
lisionai quenching (by electrons and neutral species) which removes population from the 3p^P 
state following laser excitation. Because the laser pulse is less than 5 ns long, we can measure 
this rate directly if the fluorescence lifetime is substantially longer than 5 ns. As Figure 3 
shows, the fluorescence lifetime under typical arc jet conditions is about 28 ns. The ratio of this 
lifetime to the radiative lifetime (r =35 ns) is the fluorescence quantu. . yield, about 80%. 

In summary, we have used laser induced fluorescence (LIF) to characterize the free 
stream of the 20 Megawatt Aerodynamic Heating Facility at NASA/ Ames Research Center. 

The velocities of both a major flow component (atomic oxygen) and a minor component (atomic 
copper) have been measured under a variety of operating conditions. The relative number 
density of O-atoms in the 2p^P electronic state as a function of radial position has been 
deter min ed under one set of conditions. The rate of collisional quenching for 0(3p^P) has been 
measured. In the near future we plan several extensions of this work, including absolute num- 
ber density measurements for 0(2p^P), relative number density measurements for N(2p^S), 
rotational-translational temperature measurements using 0(2p^P), and measurements of flow 
velocity as a function of radial position. A detailed comparison of the LIF results with both 
computer simulations of the arc jet flow and previous experimental measurements will be 
carried out. 


Figure 1 . LIF si gnal s for nitric oxide (in a room temperature gas cell) and atomic oxygen (in 
the arc jet) as a function of laser wavelength, with the laser at 90 degrees to the flow direction. 
(The ultraviolet excitation wavelength is half the wavelength shown, which is the dye laser 

wavelength prior to frequency -doubling.) 
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Figure 2. The same as Figure 1, but with the laser at 73 degrees to the flow direction. The 0- 
atom resonance has moved further away from the NO resonance by an amount equal to the Dop- 
pler shift. 
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normalized LIF signal 


Figure 3. Fluorescence intensity as a function of time following excitation of the 
3p3p2 i o transition in atomic oxygen with a short laser pulse. 


O-atom ftuorascanca dacay, data from Vi/33. 
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Measured and Calculated Spectral Radiation from 
a Blunt Body Shock Layer in an Arc-Jet Wind Tunnel 

Dikran S. Babikian 1 
Roger A. Craig 2 
Giuseppe Palumbo 3 
and 

Grant Palmer 4 

The shock layer spectral radiation incident on a blunt body placed in a NASA-Ames 
20A/W arc -jet wind tunnel is measured. The spectra show the atomic lines of oxygen, 
nitrogen, and copper, as well as the molecular band systems of NO , N 2 , A'+, and CN. 

The flow fields of the experiment is predicted in three separate steps. First, the equilib- 
rium flow field of the arc column and the nozzle entrance are computed from the known op- 
erating conditions of the arc-jet wind tunnel using the arc heater flow field code ARCFLO. 
Second, the nonequilibrium flow field of the nozzle and the free stream are calculated using 
the multi-temperature, one dimensional code NOZNT. Third, the nonequilibrium flow field 
of the blunt body shock layer is calculated using the axisymmetric 2d2T code, which solves 
the two dimensional Navier Stokes equations for a two temperature reacting gas Thirteen 

species; JV, O, NO, Ar, N,, 0 2 , NO*, fl+ 0* , O * . N* , Ar*, and e" are included in the 
chemistry model. 

The spectra of the incident radiation in the range 2G0nm to lOOOnm are calculated 
using the modified version of the nonequilibrium air radiation computer code NEQAIR. The 
computed and the measured spectra are compared. 

Results of analysis of the vacuum ultra violet region from 120nm to 200nm will be 
reported elsewhere. 


Experiment 


The NASA Ames 20 MW Aerodynamic Heating Facility is operated with a conical noz- 
zle. The experimental set-up, test conditions, and digital data acquisition are described in 
detail in references [lj and [2]. Figure 1 shows a schematic sketch of the nozzle and the 
spectrographic set-up. Figure 2 is a schematic sketch of the 6.0 inch flat disk model and 
the helium cooled magnesium floride MgF 2 window that transmits the surface radiative flux 
into an evacuated optical system. The incident radiation is collimated via a set of mirrors 
and focused onto the entrance slit of a McPherson spectrometer. The spectrometer has a 
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ocal .ength of 0.5 meter and operates either as a scanning monochromator or as film spec- 
trograph. The radiation from the 120nm to lOOOnm range incident on the MgF 2 window is 
measured. The data are obtained using films or photomultiplier tubes. 

Figures 3 and 4 show samples of the experimental data obtained with the photo multiplier 
tube and the film. The photomultiplier record of spectrum from 200nm to 300nm is shown in 
Fig. 3. This region is dominated by the NO-/ band system resulting from the A 2 H + — ♦ Y 2 n 
transition. Figure 4 is the densitometer trace of a film record from 500nm to 900 am Atomic 
oxygen and nitrogen emission lines are well resolved. Copper lines are seen at 324.7nm and 
327.4nm in the second and third orders. 


Computation 

" OW fc!d ° f the arC CO ‘“ m " “ d the nozzle entrance are calculated using 
tne AKCr LO (Arc Heater Flowfield) computer program [3,4,5). The inputs to the ARCFLO 

are ,, e ° P 7 ratm 5 f arameters of the 20A IW Aerothermodynamic Heating facilities given in 
Kef. 2. i he enthalpy at the entrance of the nozzle, an important parameter for solving the 
nozzle flow, ,s determined from the enthalpy distribution at the exit of the constrictor tube 
16). The nonequilibrium flow field of the nozzle and the upstream shock layer flow conditions 
are calculated using the multi-temperature nonequilibrium nozzle flow code NOZNT [71. 

igures 5 and 6 show temperature and mole fraction variation along the axis of the nozzle, 
fl c °de, 2d2T, is used to calculate the thermochemical nonequilibrium 

flowfield of the blunt body [8,9]. The full Navier-Stokes equations and 13-species chemical 
model with the latest reaction rates are incorporated. The electron excitation energy and 
he vibrational energy are assumed to be characterize by one temperature, T v , as are the 
translational and rotational energies by one temperature, T t . Figures 7 and 8 show the 
stagnation line temperature and mole fraction profiles respectively. 

The spectrum of the radiative power that strikes the stagnation point of the flat disk 

hn ITi 1S " ated USing the modified versi °n of the NEQAIR code developed by Park 
l , . j. The spectra shown in Figs. 9 and 10 are preliminary results and are obtained by 
spatially integrating the local radiative properties along the stagnation streamline. Figure 9 
p. 4 * calculated spectrum from 200nm to 300nm corresponding to measurements shown in 

ig. 3. Figure 10 shows the calculated spectrum from 500nm to 900nm that corresponding 
to measurements shown in Fig. 4. 6 

This work is an initial effort to eventually obtain a computationally based, experimentally 
va mated, description of the physical and chemical processes in an arc-jet wind tunnel. In 
the final paper, comparison of the calculations and the experimental data will be presented. 
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Figure 1. Sketch of model showing location in arc -jet wind tunnel 

flow. The stagnation region flux incident on the model aperture is 
imaged onto the spectrometer through an evacuated optical system. 




Figure 2. Sketch of the model face showing the helium cooled 

window. The view field is canted 15" from the tunnel centerline to 
exclude radiation from the arc column. 



Figure 3. Stagnation region surface flux from 200 nm to 300 nm: 
Photomultiplier output with spectrometer used in scanning 
monochromator mode. 
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Figure 4. Stagnation region surface flux from 500 nm to 900 nm: 
Densitometer trace of photographic record with spectrometer used in 
film spectrograph mode. 
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Figure 5. Calculated temperature profiles along the nozzle. 



Figure 6. Calculated mole fractions along the nozzle. 
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Figure 9. Preliminary calculated spectrum from 200nm to 300nm. 



Figure 10. Preliminary calculated spectrum from SOOOnm to tOOOOnm. 
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Abstract 

The spectral radiation from the shock layer over a 150 mm diameter flat disk model in the 
test section of an arc-jet wind tunnel is measured and computed. Three prominent molecular 
band systems, the 440-470 nm spectral range of the N 2 + First Negative (B 2 £ + u — ► X 2 E + ) ? 
the 700-780 nm apectral range of the N 2 First Positive (B 3 II 3 — *• A 3 £%j, and the 250- 9 60 
nm+spectral range of the NO-, (A 2 £ + - X 2 n) are identified and analysed. In the~case~ of 
iV 2 , transitions from u’ < 3 are observed; for NO-, band system, vibrational tranitions from 

v' < 2 are observed. These spectra are used to extract rotational, vibrational, and electronic 
temperatures. 

The experimental flow field is numerically simulated in thiee seperate steps. First the 
equilibrium flow field of the arc column and the nozzle entrance are computed from the 
known operating conditions of the arc-jet, second the nonequilibrium flow field of the nozzle 
and the free stream are computed, and third the nonequilibrium flow field of the flat disk 
shock layer is calculated. The snectra are cicuiated using the computed flow field properties 
and the nonequilibrium radiation computer code NEQAIR. The computed and the measured 
spectra are compared. 


Experiment 

The NASA Ames 20 MW Aerodynamic Heating Facility is operated with a conical nozzle. 
Two arc conditions are tested, 600 amps and 0.68 atm column pressure, and 1200 amps and 
IT atm. The spectral data are resolved spatially along the stagnation stream line yielding 
the spatial distributions of the emitting species. A McPherson spectrometer with a pho° 
tomultiplier tube is us~d to obtain the data. Our optical setup with f/50 and collection 
geometry with a de v.agnification factor of 6 limits the field of view to a cylinder with main 
axis parallel to the model face, crossing the center of the face and the jet axis. The mea- 
surement averages the radiation intensity over a cylindrical volume of 5 mm dia leter. The 
optical system forms a real image of the shock layer at the entrance slit of the n^nochrome- 
ter. The spatial profiles are obtained by moving this image over the slit. Figure 1 shows the 
NO 7 in the 250-260 nm range at 12 mm from the model surface and at the high pressure 
test condition. Vibrational tr anitions from v' < 2 are observed in contrast with previous 
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observations 1 in the freestream. Figure 2 shows the JV 2 + First Negative in the 440-470 nm 
range, and Fig. 3 shows the N 2 First Positive in the 700-780 nm range. These spectra are 
also at the same spatial location and test condition. The resulting spectral resolution is 0.2 
nm. These spectra should be adequate for extracting rotational, vibrational, and electronic 
state distribution from the approperiate species and comparing with the computations. 

Computation 

The equilibrium flow field of the arc column and the nozzle entrance are calculated using 
the ARCFLO (Arc Heater Flowfield) computer program. 2 ’ 3 ' 4 The inputs to the ARCFLO 
are the operating parameters of the 20 MW Aerothermodynamic Heating facilities. The en- 
thalpy at the entrance of the nozzle, an important parameter for solving the nozzle flow, is 
determined from ‘.he enthalpy distribution at the exit of the constrictor tube . The nonequi- 
libiium flow field of the nozzle and the upstream shock layer flow conditions are calculated 
using the multi- temperature nonequiiibrium nozzle flow* cooe NOZNT. J 

The two-dimensional code, 2d2T, is used to c; 1-ulate the thermochemical nonequilibrium 
flowfield of the model 6 ’ 7 . The full Navier-Stokes equations and 13-species chemical model 
with the latest reaction rates are incorporated in this code. The *-'.ectron excitation energy 
and the vibrational energy are assumed to be characterized by one temperature, T v . Similarly, 
translational and rotational energies are characterized by one temperature, T) . Figures 4 and 
5 show a typical such solution for the stagnation line temperature and mole fraction profiles 
respectively. 

The spectrum of the emitted radiative power at a given spatial location away from the 
surface of the model is calculated using the modified version of the NEQAIR code developed 
by Park 8 ’ 9 using the computed properties. 

This work is an initial effort to eventually obtain a computationally based, experimentally 
validated, description of the physical and chemical processes in an arc-jet wind tunnel. In 
the final paper, comparison of the calculations and the experimental data will be presented. 
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Abstract 

Radiation from the nitric oxide band systems emit- 
ted by the flow in the test section of a 20 MW arc-jet 
wind tunnel was measured and computed. The settling 
chamber pressure and enthalpy were 2.4 atm and 28il0 
MJ/kg, respectively. The measurements were made with 
a 0 3 meter McPherson spectrograph using photographic 
films in the wavelength region from 225 nm to 305 nm 
Of the four band systems of NO (/?, 7, 6, and r), vi- 
brational transitions were observed from only the upper 
v' = 0 levels. Excitation temperatures were deduced by 
comparing the expert ntal spectrum with those calcu- 
lated using the nonequ. urn radiation code NEQAIR. 
The rotational, vibratio 1, and electronic excitation tem- 
peratures deduced from the data were: T r = 560 ± 50° 
K, T v < 950 ±50° K, and T, x = 11500 ±520° K, respec- 
tively. A multi- temperature nonequilibrium nozzle flow 
code NOZNT was used to calculate the nozzle flow. The 
calculated temperatures were T r = 560° K, T„ = 950° 
K, and electron thermal temperature T e = 6100° K, re- 
spectively at 30 MJ/kg. The sensitivity of the calculated 
temperatures to the uncertainty in the measured arc-jet 
enthalpy is examined. 

Introduction 

For the last three decades arc-jet wind tunnels have 
been used to simulate the high enthalpy environments 
encountered during hypervelocity atmospheric entry of 
space vehicles. The heat shields of all such vehicles, in- 
cluding those for the Apollo and the Space Shuttle vehicle, 
have been develope d by testing in arc-jet wind tunnels 
Testing of heat shield materials will become even more 
critical for the proposed Mars mission because the entry 
velocity will be higher than those for the past missions. 
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In spite of the importance of the arc-jet facilities, the 
aerothermodynamic states of the flows produced in them 
are not yet well understood. This is due to the complex- 
ity of the physical and chemical phenomena occurring in 
these facilities. The flow undergoes vibrational excitation, 
dissociation, and ionization. Because of the relatively low 
density environment produced in'fche facility, it is gener- 
ally in nonequilibrium in both chemical composition and 
internal modes. 

There have been some efforts to characterize the arc- 
jet flows by using the available experimental and com- 
putational techniques. 1,2 Despite such efforts, many un- 
certainties still exist. In particular there has not been a 
measurement of vibrational or t -tronic temperatures of 
nitric oxide (NO). These temperatures of NO are impor- 
tant because NO plays an important role in the chemistry 
of the nozzle flows. 3 

The objectives of the present work are to experimen- 
tally characterize the behavior of NO and to explore 
whether it can be numerically reproduced. Radiation 
from NO was measured in the test section stream of a 20 
MW arc-jet wind tunnel at NASA Ames Research Center. 
From the data, we determined its rotational, vibrational, 
and electronic excitation temperatures that are consistent 
with the limited data. The results of the measurements 
were compared with the theoretical calculations. These 
calculations included the computation of arc constrictor 
flow, nozzle flow, and nonequilibrium radiation. Due to 
the complex nature of the flow in the settling chamber, 
accurate prediction of the flow properties in this chamber 
remains uncertain and is under investigation 4 The calcu- 
lated temperatures and the synthetic radiation spectrum 
are generally in fair agreement with the measurements. 

Experiment 

Test Conditions and Experimental Setup 

The arc heater of the 20 MW arc-jet wind tunnel at 
NASA Ames Research Center used in the present work 
(Aerodynamic Facility) consists of a constrictor tube that 
is 216 cm Ion J with a 6 cm internal diameter. The con- 
strictor is followed by a settling chamber 10.4 :m in diam- 
eter and 25.6 cm in length. From the settling chamber, 
the flow is expanded through a convergent-divergent noz- 
zle of throat diameter 3 8 1 cm, to an exit diameter of 
45.72 cm with an 8° half angle and an area ratio of 130. 

The arc heater was operated with a settling chamber 
pressure of 2.4 atm, a flow rate of 0.0864 kg/sec, electric 
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current of 2000 A and a voltage of 2600 V. The test gets 
consisted of 87% air and 13% argon by mass. The average 
value of the flow enthalpy based on heat balance measure- 
ments. i.e. electrical energy input minus the energy lost 
in the cooling water divided by the mass flow rate, was 22 
MJ/kg. The {low enthalpy at the nozzle centerline based 
on heat transfer rate measurements to a spherical blunt 
budy 5 was 28 MJ/kg. 


ADJUSTABLE 
CONCAVE MIRROR 




Fig. i. Schematic drawing of the experimental set-up. 

Spectral measurements were made in the freestream 
flow 35.6 cm downstream of the nozzle exit. Because of 
the low density nature of the expanding flow, the radia- 
tion emitted by the gas was very weak. Therefore in or- 
der to obtain sufficient signal, a McPherson spectrograph 
model 218 with a small f-nurnber and a short focal length 
was used to collect the data. The spectrograph has an 
f-number of 5,3 and a focal length of 0.3 meters. This 
spectrograph is known as a criss-crossed Czerny-Turner 
optical system 6 . It is an unconventional mounting ar- 
rangement where the light beam from the two mirrors 
and the plane grating are crossed in order to obtain a 
low f-number. The signal was recorded on a film on this 
instrument by an exposure of 30 sec. A Joyce Loebel mi- 
crodensitometer model MK III CS was used to measure 
the optical densities recorded on a High Speed Infrared 
Film (HIE). The microdensitometer traces were digitized 
and transferred to a computer to be calibrated, analyzed, 
and compared with the theoretical predictions. 

Further description of the experimental setup, instru- 
mentation, data collection, and data calibration, are 
presented in References 7 and 8. Figure 1 shows the 
schematic sketch of the experimental optical set-ui 


Experimental Spectra 

As can be seen from Fig. 1, the radiation observed by 
the spectrograph is an integration of the intensity along 
the optical path across the flowfield. How'ever since radia- 
tion power is generally strongly dependent on the control- 
ling temperatures, and since these temperatures are the 
highest along the centerline of the nozzle flow, it is likely 
that the observed spectra are representative of those in 
the centeriine region of the nozzle flow 

The densitometer trace of the record taken in the spec- 
tral range 450 nm to 850 nm, in the first order is shown 
in Fig. 2. According to the spectrum shown in Fig. 2, 
the only significant radiation emanates from NO7 and 6 
band systems (in second order). 



Fig. 2. Densitometer trace of the 450 nm to 850 nm (first 
order) spectral range. 

The measured intensity was corrected for the film re- 
sponse in two steps using a National Bureau of Standards- 
certified calibrated lamp. First the film was calibrated for 
the relative intensity at specific wavelengths using neutral 
density step- wedge filters. Second the film was calibrated 
for the film sensitivity, the spectrograph grating, and the 
entire optical set-up against wavelength. Figure 3 shows 
the corrected intensity trace of the NO band systems ad- 
justed to the first order in the 225 nm to 305 nm wave- 
length range. 

As indicated in Fig. 3, the spectral bands detected 
by the spectrograph were identified as the NO7 and NO<§ 
band systems in second order. The NO7 band system 
corresponds to the A 2 £ + — ► X 2 II transition where A 2 E + 
is an upper electronic state and X 2 II is the ground elec- 
tronic state of the NO molecule. The NO<5 band system is 
produced by transitions from the C 2 II — X 2 H electronic 
states. The NO 0 band system (B 2 II+ — X 2 II), which is 
in the same spectral region as the NO7 band system (see 


AIAA 93-2800 3 


Fig, 4), was not present in the recorded spectrum. For 
both the 7 and 6 systems only the vibrational transitions 
from the upper, v' = 0 , level were present: i.e. ( 0 , 1 ), (0 2 ) 
(0,3), etc. V 



Fig. 3. Corrected intensity trace of the 225 nm to 305 nm 
spectral range. 



INTERNL CLEAR CiSTANCc -nm) 

Fig. Potential curves of the NO molecule. 


Temperature Determination 


The excitation temperatures of the flowfield were de- 
duced by comparing the relative intensities in the recorded 
spectra with those calculated using the latest version of 
the nonequilibrium radiation code NEQAIR. 9 10 The code 
was used to calculate the emission spectra in the 225 nm 
to 305 nm range. The thermochemical model used in thi 
code has been well documented in Refs. 11-14 and there- 
fore will not be described here. The nonequilibrium ex- 
citation portion of the code was bypassed by specifying 
Boltzmann distributions of electronic, vibrational, and ro- 
tational excitation modes. Also, the radiation was very 
weak and the gas was assumed to be optically thin. 

The width of a vibrational band is dictated mostly by 
the rotational temperature of the molecule. In the present 
work, the rotational temperature was determined from the 
(0.2) band of the 7 system. The temperature was altered 
in the calculation until the calculated value of the band 
width at half its maximum intensity matched that of the 
experimental spectrum. The rotational temperature so 
determined was T r = 560 ± 50° K. The source of the un- 
certainty ±50° K is attributed to the uncertainties about 
film characteristics . 1 8 The ratio of the relative intensity of 
the 6 system to that of the 7 system depends on the elec- 
tronic excitation temperature governing the two systems. 
Through a computational process similar to that used for 
the rotational temperature, an electronic excitation tem- 
perature of T er = 11500 ±520° K was obtained. In Refs. 

7 and 8 , T ei was erroneously stated to be 7560±340° K. 
This error was caused by the use of the old values of inten- 
sity parameters (transition moments). To determine the 
vibrational temperature of NO, radiation from two vibra- 
tional levels is needed. However, only bands from v' — 0 
state were observed in the data. From the fact that the 
radiation from the 1 / = 1 state of the upper electronic 
state of NO 7 was absent in the observed spectrum, only 
the upper limit of the vibrational temperature of NO can 
be determined. In order for the radiation from the 1 / = 1 
state be equal to or less than the observed noise level, 
the vibrational temperature was determined to be bellow 
950±50° K. 


Figures 5a and 5b show the calculated spectra of the 
radiation emitted by the NO band systems. In Fig 5 a 
the emission spectrum from the 7 and S systems is shown. 
This spectrum is in good agreement with the experimen- 
tal spectrum shown in Fig. 3. The emission spectrum 
from the 0 and e systems is shown in Fig. 5 b. Note that 
the radiation from these systems is an order of magnitude 
weaker than that from 7 and 6. Also, the /3 and c systems 
are not present in the recorded spectrum shown in Fig. 3 . 
The absence of the 0 band system indicates that the B 2 T1 




state is not populated o v .. 

though it is located at nearly the same energy level as the 
state from which the 7 system emanates (see Fig. 
4). The B 2 II state has a larger internuclear distance than 
the A 2 D+ state. The electronic states of the NO molecule 
may be populated through a selective populating process 
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according to internuclear distance under nonequilibrium. 

The upper state of the NOc band system (D 2 S + ) lies 
very close to the the upper state of the NO7 band system 
as shown in the energy potential curves in Fig. 4. How- 
ever, the strong bands of N Oe lie in the vacuum ultra- 
violet and only a few of the weaker bands are in the spec- 
tral region of our experiment. These weak vibrational 
bands (0,5), (0,6), and (0.7) lie under the much stronger 
vibrational bands ( 0 , 0 ), ( 0 , 1 ), and ( 0 , 2 ) of the NO 7 sys- 
tem and their contribution to the measured radiation is 
insignificant. To verify the existence of NOe in the free 
stream flow, further emission measurements are needed in 
the vacuum ultra-violet spectral region where NO7 band 
system does not exist. 
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Fig. 5. Calculated spectrum of NO; a) 8 and 7 band 



Flowfield Calculations 
Nozzle Inlet Conditions 

In order to carry out the nozzle flow calculations, one 
must know the conditions at the nozzle entrance. The 
entrance condition is dictated by the flow processes in 
the settling chamber Since the settling chamber receives 
its flow from the arc constrictor, one must first calculate 
the constrictor flowfield. The flow conditions at the exit 
of the constrictor were calculated using the arc heated 
flowfield code ARCFLO. 16,1 ' The code assumes the flow 
to be in thermochemical equilibrium and uses real gas 
properties . 16,18 





Fig. 6. Radial distribution of velocity, temperature, and 
enthalpy profiles. 

Figure 6 shows radial distributions of velocity, tem- 
perature, and enthalpy at the exit of the constrictor. The 
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centerline values were 225 m/sec, 12200° K, and 60 MJ/kg 
respectively. The calculated mass averaged total enthalpy 
at the exit of the constrictor was 22 M J /kg. The calcu- 
lated averaged enthalpy value agrees with the value deter- 
mined from heat balance. However, the calculated center- 
line enthalpy value is much greater than the value based 
on heat transfer rate. 

The discrepancy between the two centerline enthalpy 
values is believed to be due to the circulation and mix- 
ing between the hot gas in the centerline region and the 
cold gas in the peripheral region. The resulting distribu- 
tion of these properties in the settling chamber is diffi- 
cult to predict, although an effort is presently being made 
to analyze it. 4 However, because of this phenomenon, it 
is expected that the distribution of temperature and en- 
thalpy at ’ts exit will be more uniform than that shown in 
Fig. 6. Therefore it is reasonable to assume the centerline 
enthalpy to be in the range of 28dbl0 MJ/kg 

In the present work, calculations were carried out for 
enthalpy values between 20 and 60 MJ/kg in an incre- 
ment of 5 MJ/kg. The value of 30 MJ/kg, which is close 
to the 28 MJ/kg value deduced Lorn heat transfer mea- 
surement, was used as the inlet enthalpy for the nozz-e 
flow calculation. 




Fig. 7. Calculated velocity and pressure variation along 
the nozzle axis. 

Nozzle Flowfield 

The computer code NOZNT 2 15 was used to calculate 
the flow along the centerline of the nozzle of the arc-jet 


wind tunnel. The code solves a one-dimensional steady 
flow through a convergent-divergent nozzle in the dissoci- 
ated and ionized regime. In the nozzh flow calculations, 
the nozzle geometry was assumed to be hyperbolic, that 
is, conical with a smooth transition at the throat, of the 
form A/A* = 1+cx 2 It assumes thr entrance and the be- 
ginning section of th nozzle to be at equilibrium. The 
rest of the nozzle is solved assuming multi-temperature 
nonequilibrium flow. The equilibrium portion of the flow 
is solved with the pressure-specified method. 19 

The nonequilibrium calculations of the flow are started 
in the converging portion of the nozzle upstream of the 
throat. The nonequilibrium calculation also uses the 
pressure-specified method up to the point where the 
frozen Mach number is 1.8. From this point the area- 
specified method is used. 2 In this region of the flow the 
equations of conservation of species, vibrational energy, 
and electron-electronic energy are solved numerically. 

Figure 7 shows the calculated pressure and velocity 
along the axis of the nozzle. The code assumes that the 
translational temperature T and rotational temperature 
T r are the same. Likewise, the electron translation tem- 
perature T« and electronic excitation temperature T cr are 
taken to be the same. However, molecular species i are 
allowed to have different vibrational temperatures T Vi . 


a> 

3 

aj 

v_ 

o 

Cl 

E 

© 


1 


1 




I 



n 1 — i — T — i — 1 — r 

0.0 0.4 0.8 1.2 


Distance from entrance, m 


Fig. 8. Calculated excitation temperatures along the noz- 
zle axis. 


The reactions and excitation rates used in the code 
and solution procedure are discussed in Ref. 15. Those 
rate parameters were used in the present work except for 
those for ionic recombination of N + and 0 + . Numerical 
difficulties were encountered for some calculations where 
enthalpy was higher than 45 MJ/kg. These difficulties dis- 
appeared when the ionic recombination rates were multi- 
plied by 10. For the cases where there were no numerical 
difficulties, the change of the rate constants resulted in 
negligible differences. Since these rate constants are un- 
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certain by about an order of magnitude the changes are 
acceptable for the purpose of the present work 



Distance from entrance, m 
Fig. 9. Calculated mole fractions along the nozzle axis 

Molecular oxygen was excluded from the calculations 
because it is almost entirely dissociated at the test con- 
ditions and its concentration is extremely small. Figures 
8 and 9 show the temperature and mole fraction profiles 
along the nozzle axis. At the point where the ex. 'rimen- 
tal measurements were made, the calculated temperatures 
were T r = 560° K, T„(NO)=950® K, and T e = 6100° K 
at the enthalpy value of 3U MJ/kg. 

Discussion 

The excitation temperatures deduced from the experi- 
mental spectrum are shown in Fig. 8 by the circle symbols 
with error bars and are compared with the NOZNT calcu- 
lation. The figure shows that the calculated translational- 
rotational and the NO vibrational temperatures are in 
good agreement with the corresponding experimentally 
deduced values. The difference between the measured 
« and the calculated T, may be due to the breakdown 
of the assumption that the two temperatures are equal 
However both these temperatures. T, and T„, are much 
igher tnan T t .(NO). and thereby it underscores the need 
or a multi-temperature description of the flow. The vi- 
brational temperature of N 2 was not measured in this ex- 
periment However, the accuracy in the calculated value 
on„(N 2 ) was verified in Ref. 2 by comparing it with 
experiments where N 2 vibrational temperature was mea 
sured. 

In order to examine the sensitivity of calculated tem- 
peiatures to the uncertainty in the measured arc-jet en- 
t a py, the foregoing calculations were repeated for differ- 
ent values of nozzle entrance enthalpies. The results a'e 
shown m Figure 10. It can be seen from the figure that, as 
the enthalpy increased from 20 MJ/kg to 60 MJ/kg, the 


translational-rotational temperature increased from 550® 
K to 1500° K. However, the electron temperature reached 
a maximum value of 6200® K at H=35 MJ/kg and de- 
creased thereafter as the enthalpy increased further. At 
enthalpies greater than 45 MJ/kg the electronic tempera- 
ture and all vibrational temperatures are in equilibrium 
whereas at enthalpies below 25 MJ/kg the translational 
temperature and the NO vibrational temperature are in 
equilibrium. 





Fig. 10. Effect of inlet enthalpy on nozzle flowfield. 

The value of T„(NO) found in the present experiment 
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indicates that the the enthalpy is below 35 MJ/kg. A pre- 
vious work 20 indicates that the enthalpy value deduced 
from the heat transfer measurement tends to underesti- 
mate the enthalpy in typical arc-jet conditions. Therefore, 
the enthalpy is most likely between 28 MJ/kg, the value 
determined from the heat transfer measurement, and 35 
MJ/kg. It is to be noted that at the enthalpy of 60 MJ/kg, 
which is the value ot the enthalpy at the constrictor cen- 
terline, the concentration of NO is very small and it is 
unlikely that its radiation could be observed experimen- 
tally. 

As shown in Ref. cO. enthalpy can be best determined 
spectroscopically by observing radiation emitted from the 
shock layer formed in front of a blunt body. Such effort 
is currently in progress. Measurements of the spectral ra- 
diation (120 nm to 900 nm) incident on the stagnation 
region surface of a blunt body in the test section have 
been made 4 **’- 2 and the analysis described in this paper is 
being extended to include the shock layer modeling. Pre 
liminary predictions of the incident radiation spectrum 
have been made and are being compared with the exper- 
imental data. 20 * 24 Results of this analysis will be used to 
define future arc-jet tests to continue development of arc- 
jet wind tunnel modeling. A validated model of these fa- 
cilities will help design advanced arc-jet and extend their 
usefulness as aerothermodynamic testing facilities. 

Conclusions 

The presented results show that by using the centerline 
enthalpy value deduced from heat transfer measurement 
and the NOZNT code, one can predict the free stream 
conditions in an arc-jet wind tunnel flow fairly well. The 
translational-rotationai temperature and the vibrational 
temperature of NO can be closely reproduced by NOZNT. 
The calculated electron-electronic temperature T e is ap- 
preciably lower than the measured electronic excitation 
temperature T er of NO. Compared with the electron and 
electronic temperatures, the vibrational temperature of 
N 2 and NO are significantly lower at enthalpies less than 
45 MJ/kg. In the same enthalpy range, the vibrational 
temperature of NO is significantly lower than that of No 
The ii„halpy deduced from the spectroscopic measure- 
ments agrees approximately with that deduced from heat 
transfer measurement. 
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Abstract 

Current and future calculations of nonequilibrium 
shock layers require the use of a very large number of 
equations, due to a multiplicity of chemical species, 
excited states, and internal energy modes. The com- 
putational cost associated with the use of standard 
implicit methods becomes prohibitive; it is, there- 
fore, desirable to examine the potential of several 
methods and determine if any can be projected to 
be more efficient and accurate for large systems of 
equations. In this paper we examine the performance 
of several implicit schemes on som* simple practical 
examples of reacting flows. The Euler equations are 
solved by three different implicit methods, and two 
methods of coupling between the fluid dynamics and 
the chemistry are studied. Several cases of stiffness 
are considered, and both one and two-dimensional 
examples are computed. We conclude on with some 
remarks on the accuracy, stability and efficiency of 
these various methods. 

I. Introduction 

The modern Computational Fluid Dynamics 
(CFD) tools are becoming increasingly useful in com- 
puting complex flow conditions, which generally in- 
clude non-equilibrium phenomena. There is a gen- 
eral need for increasingly complex modeling of the 
thermo- chemical properties of the gas, and for the 
modeling of larger systems. For example, the model- 
ing of shock layers around ablating bodies requires a 
very large set of chemical species and chemical reac- 
tions. Although some approximate formulations can 
be used in the preliminary design phase of space ve- 
hicles or experiments, the modeling of the complete 
kinetics is desirable or even required when the non- 
equilibrium effects become dominant: this happens 
for example as the flow expands around the shoul- 
der of a vehicle, or when the object is reduced in 
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size for insertion into an experimental facility. The 
situation can be further complicated due to the fact 
that most flows realized in ground-based experimen- 
tal facilities are themselves not in chemical or ther- 
mal equilibrium. Another example concerns highly 
ionized and radiating flows, which are likely to be 
found at high re-entry velocities, or their equivalent. 
It may be required, for these cases of plasma condi- 
tions, to account for non- Boltzmann distribution of 
the excited states. This problem may also require 
us to model the plasma with a complete collisional- 
radiative model of the plasma, and convect all the 
excited states, forcing us to use a large number of 
equivalent species. In addition, the internal relax- 
ation processes will be locally stiff and difficult to 
model. These upcoming challenges in CFD technol- 
ogy will require the development of efficient meth- 
ods for a very large number of species, and for pos- 
sibly stiff couplings to complex internal processes. 
Since we want a method that allows us to reach 
the steady state with minimal computational effort, 
it seems desirable to jse an implicit method. On 
the other hand, since we may need to couple the 
hydrodynamics to several other physical processes 
(collisional-radiative processes, radiation transport, 
electro-magnetic couplings, etc..), we may favor the 
use of the Operator-Splitting (OS) method. The lat- 
ter must be contrasted with the Fully-Coupled (FC) 
approach, which attempts to provide a more accu- 
rate and more stable way to couple the different pro- 
cesses. It is not clear a priori which method is more 
accurate, stable, efficient, or practical; the use of one 
method versus the other may depend on the type of 
flows beir f computed, the type of computer archi- 
tecture used, or even the personal preferences of the 
modeler. During the course of the present work, we 
will compare the FC and OS methods (for chemistry) 
and therefore add fuel to the debate. Our search for 
an efficient numerical method, extendable to large 
systems, will also include the effect of required grid 
accuracy on the solution, and its impact on the effT 
ciency of the numerical approaches used. 


II. Numerical Methods 

II- A: The N x N Biock-Tridiagonal method 

The Euler equations describe the convective pro- 
cess. and are written (in ID) as: 
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where E , H are respectively the total energy and to- 
tal enthalpy, per unit volume. The internal energy 
is Ei = P* f T dT'CvsiT'), and the gas mixture 
follows the ideal equation of state: 


implicit system for 6 = 1, while second-order time 
accuracy is obtained for 9 - 1/2. Second-order spa- 
tial accuracy is achieved by evaluating the fluxes at 
the cell interfaces i ± 

F,±i = ^(F, + F,±i) 

A final modification to the fluxes is made to assure 
monotonicity. The Euler system is sr hyperbolic sys- 
tem, and has a set ci’ real eigenvalues (characteristic 
speeds). The Jacobian can then be written in the 
form: 

A = T _1 AT (4) 

where the matrix of eigenvalues 
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P = NkT =(y- 1 )E, = (7 - l)(E - ipu 2 ) 


The individual species densities are denoted by p,, 
and p = P* ■ This formulation is for a single fluid 
(one mass-averaged velocity), in thermal equilibrium 
(one temperature). This equation is discretized over 
a finite size mesh to yield the form: 
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The subscripts i, i ± \ indicate that the variables 
are evaluated at computational cells (center) i, and 
cell interfaces ii 1. V and S are respectively the 
cell volumes and surfaces. This finite- volume for- 
mulation will be used throughout the paper. The 
equation has been discretized in time as well, and 
the expression AQ describes the difference between 
the flow variable evaluated at two time levels, n and 
n+ 1. The expression on the RHS of equation (2) 
must be further specified: the fluxes are a function of 
Q, and can be evaluated at time level n + 9 through 
the linearization approximation: 

~ +0A(Q)(AQ) i+i (3) 

where A = is the Jacobian matrix. The explicit 
Euler system of equations is obtained for 0 = 0, the 
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is diagonal and real only, T,T _1 are transfer matri- 
ces between the space of primitive variables Q and 
‘characteristic’ variables, and c is the speed of sound. 
The spectrum of eigenvalues can be split into and 
positive and negative values, which indicate the di- 
rection of flow of the characteristic variables at the 
cell interface. The flux at an interface can now be 
written as: 


F (n + ») ^ F (^ +9A + ( Q )A Q. +5A - ( Q )AQi + 1 (6) 


where 

A* =T _1 A* T 

and A* is the set of eigenvalues which are respec- 
tively positive (negative), zero otherwise. Using 
this formulation, the discretized version of the Eu- 
ler equations becomes: 


MAiA+^AQ.-i 

+[1 + 0A*A+ +i - flAfA- ^AQ. (7) 

The RHS of equation (7) can be modifed for mono- 
tonicity, while conserving its second-order accuracy 
in space. The technique used throughout this work 
follows closely the TVD method of Harten [1]. 
Greater stability is generally obtained if the implicit 
LHS of equation (7) has its spatial accuracy reduced 
to first order. This consists in evaluating the Jaco- 
bian matrices at the cell centers, according to the 
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characteristic flow directipn. The final version of the 
system is obtained for the fully implicit case (0 = 1); 

[-A+ l A<]AQ J ., 

+[1 + A, + A<- A-A<]AQ, (8) 

+[A- +1 A<]AQ I+I = RH& n) 

This is a tridiagonai system of N x A r block ma- 
trices, where N = N $ 4* D + 1, N , is the number 
of species, D is the spatial dimensionality. Solv- 
ing this system can be done by the standard tech- 
nique of gaussian elimination and back substitution, 
with LI decomposition of the block matrices (see 
for example [2]). This requires that the matrices 
that compose the diagonal band be inverted twice 
for each grid point. It turns out that the cor- 
responding algorithm has a number of operations 
that scales as N 3 . This technique is called the 
N x N Block-Tridiagonal solve*, and will be used in 
this paper. 

A similar approach can be used for 2-dimensional 
flow, leading to a pentadiagonal system of N x N 
bock matrices. Another approach, which consumes 
less memory and has a lower number of operations, 
uses the technique of dimensional splitting [3]: we 
effectively solve the tridiagonal system twice, once 
for each direction^. This is the approach used here. 

I I-B: The Scalar Tridiagonal method 

It is clear that as the number of species grows, 
the N 2 dependence of the computational load will 
rapidly make this technique impractical. It is de- 
sirable then to search for a method that scales less 
rapidly with the number of species. One such method 
can be easily obtained by considerable simplification 
of the matrix structure. Note that the split Jacobians 
are bounded respectively from below and above: 

A + = T -1 .A + .T < max{A} 1 (9) 

A~ = T -1 .A“.T > min{A] • 1 (10) 

where we have used the (signed) maximum and min- 
imum eigenvalue present in A. The Jacobian matri- 
ces, thus replaced into equation (8) are proportional 
to the unit matrix. There is only a scalar opera- 
tion to perform, instead of a full block matrix in- 
version. The number of operations now scales as N: 
this scheme is called the Scalar-Tridiagonal solver, 
and will be compared to the previous one. 


^This splitting is a form of Operator- Splitting. 


There are several disadvantages to the scalar 
technique: the first, loss of time accuracy, is not of 
immediate relevance, since we are mostly concerned 
with the achievement of steady-state. The second 
is a loss of accuracy: this is specially of concern in 
subsonic regions, where the spectrum of eigenvalues 
is originally very different from the maximum ( nin- 
imum) value. In supersonic or hypersonic regions, 
this is not a problem, since u >> c, and A cr ul . i.e 
the spectrum is nearly scalar. We may expect there- 
fore some loss of accuracy, or even stability, when 
using the scalar method. 


II-C: The iVj-Split Tridiagonal method 


We will also investigate another method, based on 
the formulation of a multi-fluid system of equations. 
Let us consider the following system of equations: 

* U : ) + £ ( * 5:“’ ) =0 (11) 


and similar systems for other species. In this for- 
mulation, each species is attributed its own momen- 
tum density and energy density. In the limit of very 
strong coupling between the momentum and energy 
densities of each fluid component, we can enforce 
a unique velocity and unique temperature for this 
multi-fluid description. If we were to solve each sub- 
system by the block tridiagonal method, we would 
req lire the inversion of a 3 x 3 (4 x 4 in 2 dimensions) 
matrix. We repeat the method for each species, and 
the overall cost now scales as N , . For a large num- 
ber of species, we expect a cross-over between this 
method and the N x N block method, by compar- 
ing the costs; for example, 4 3 N $ versus (N g + 3) 3 . 
However, a further reduction in cost can be achieved 
with the following approximation. Assuming that 
all species have nearly equal molecular masses, and 
that their individual specific heats are nearly equal, 
we can replace the partial pressure: 

p, = n,RT ~ —P 
P 


and use a constant average 7 in the formulation of 
the derivatives which compose the Jacobians. These 
Jacobian matrices become then identical. This has 
a rather drastic effect: the block matrix inversions 
need to be performed once, instead of once for each 
species. This lowers considerably the overall CPU 
requirement. This formulation of the solver will be 
called Aj-Split Tridiagonal solver, and will compared 
with the tw' previous 
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II-D: Hydrodynamics-Chemistry Coupling 

In the equations considered so far, the source 
terms on the RHS are non existent; for a reacting 
gas, there will be a source term W which operates 
on the species densities only. The total energy now 
includes the energies of formation: 



and is not not affected by the chemical reactions 
the internal energy is obtained from the conserved 
total energy, after subtracting the kinetic and for- 
mation energies. The change in formation energy of 
the mixture, due to chemical reactions, will thus be 
converted into a change of temperature. 

The chemical source term can be treated im- 
piicitely as weil: if H is the corresponding Jacobian 
for the source term (fl = ^), solving for the chem- 
istry alone would read as: 

[l + n i At]A ch Q = W, (n) A« (121 

By solving separately for the convective and chemical 
terms, one obtains two variations at the end of the 
computational step: the global variation will then 
be a direct sum of the the contribution from each 
process. 

AQ = AcvQ + A C hQ (13) 

This procedure is called* the Explicit-Coupling 
method (EC). 

We see that in effect, we have split the com- 
putation in two parts, for each physical process. 
For that reason, this method is also called the 
Operator-Splitting (OS) method (see for example 
[3]). Another form of operator splitting consists in 
using the change induced by one process as a starting 
point for the other process: a temporary solution Q 
is used, such that: 

Q = Q< n) + A cv Q (14) 

Q( n +l) = Q + AchQ (15) 

where now the change induced by chemistry is ob- 
tained by using the modified solution Q in the ex- 
pression of the source term and Jacobian, W,fi in 
eq. (12). This formulation of the Operator-Splitting 
method is based on fractional steps, and is best de- 
scribed in [4]. Since we will examine both meth- 
ods, we will reserve the term Explicit-Coupling (EC) 

t Also called the Loosely- Coupled approach. 


for the method described in eq. (13), and the term 
Operator-Splitting (OS) when using the method de- 
scribed by eqs (14) and (15). 

Another approach is to solve for both the con- 
vective and chemical processes simultaneously The 
chemical Jacobians can be brought into the LHS, and 
equation (8) is modified to: 

[“A+.jAfjAQ,-! 

+[1 + A+ At - A," At + f^AtjAQ, (16) 
+[A- +1 A<]AQ,+i = RH 5< n) 

where now the RHS includes the evaluation of the 
chemical source terms at time level (n). This method 
is called the Implicit, or Fully-Coupled (FC) ap- 
proach. The H matrix is dense, and the FC approach 
described above is possible only when combined with 
the N x N Block-Tridiagonal solver. Including it in 
the N t - Split solver would require serious modifica- 
tions, and has not been attempted here. Similarly, 
by approximating the Jacobian H by a scalar (using 
again the maximum eigen va ue), one could use the 
FC approach with the Scalar Tridiagonal solver. We 
found that in many cases this approximation usually 
leads to very poor results for the chemistry, and will 
not be used here. 

II-E: Chemistry Sub-Cycling 

There are additional modifications one can make 
when using the OS or EC approaches: since the fluid 
dynamics and chemistry are computed separately for 
a global time step A f, one has considerable flexi- 
bility in the methods used for each process. No- 
tably, the accuracy of the chemistry can potentially 
be improved by sub-iterating (more precisely sub- 
cycling) the chemistry by using smaller time steps 
it. This may be required to improve the accu- 
racy, because the chemical reactions are non-linear 
processes: linearization errors become important in 
some highly non-equilibrium situations. The cou- 
pling of the chemistry to the temperature can also be 
estimated at each sub-step, by looking at the change 
induced in the average formation energy of the mix- 
ture. When the chemistry is sub-iterated (SI), the 
global variation is obtained by using eq. (13), but 
when the global change due to chemistry is obtained 
as follows, using sub-iterations (m = 1,2, . . .) of the 
chemistry: 

AchQ (m + i) _ A ch Q (m ' + [1 + nit]- 1 W 6t (17) 

Finally, the coupling of the chemistry to the con- 
vection can also be computed at each sub-step. For 
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example, after computing the variation A CV Q using 
one of the Tridiagonal solver listed above, one can lin- 
earize it during the global time interval. The global 
variation (for both processes) is then obtained as fol- 
lows: 

AQ {m+l ' = AQ (ml + [1 + fl<5<] _1 W6t (18) 


This formulation replaces equation (13). We will 
reserve the term Sub-Iterated Coupling (SIC) to de- 
scribe this particular form of the splitting method 
between chemistry and convection. 

II-F: Performance 

The methods used can be classified, according 
to the treatment of the convective process and the 
method of coupling with the chemistry: the designa- 
tion of the methods which will be studied are listed 
in Table 1. 

The relative performance of all schemes is demon- 
strated in Figure 1. Figure 1-a (top) shows the CPU 
spent (per iteration and per grid point) by an im- 
plicit method, normalized by the same quantity for 
the explicit method. Since the explicit method scales 
almost exactly as the number of equations, both the 
Scalar and the JV,-Split method will show a nearly 
flat behavior when plotted versus the number of 
species. This is confirmed in Figure 1-a. Notice also 
that the relative cost of the Scalar method is very 
sm<dl, while the N x N Block method climbs very 
rapidly: the latter is still quite expensive, even for a 
small number of species. The leftmost data point at 
N, = 5, for example, shows that the N x N Block 
method is 10 times more expensive than the explicit 
method. Although this number is not an absolute, 
and can be reduced after a strong effort in code writ- 
ing (by ‘hard-wiring’ the operations, for example). 
At best, this time may be reduced by a factor of two. 
Still, the conclusion is inevitable: as the number of 
species grows, the implicit scheme is efficient only if 
it can be operated at large CFL numbers. Practically 
speaking, stability limitations will limit the CFL to 
the neighborhood of 4-5. These stability problems 
arise from transient phenomena, dimensional split 
errors and/or coupling errors with the chemistry or 
other internal relaxation processes, iiigher values of 
CFL number can potentially be achieved when the 
flow is very close to the steady state and when the 
flow is non-stiff. Since we are mostly interested into 
reaching the steady state (and having to go through 
the transients) and into stiff problems, this is of lit- 


tle interest to us. These limitations will be demon- 
strated on some practical sample cases in the next 
sections. 

III. One-Dimensional Shock 

As a first test case, we will model the propagation 
of a 1 -dimensional shock, from an impulsive start 
This case will mimic the establishement of a two- 
dimensional sho^k layer around a blunt body. WV 
use a grid of 200 cells, evenly spaced, with a per- 
fectly reflecting wail on the right hand side. The 
flow is incoming from the left at high velocity, and 
impinges on the wall. A shock is created at the re- 
flection and propagates back upstream into the hy- 
personic flow. Although strictly speaking this flow 
is unsteady, the profiles become steady in a frame 
attached to the shock. The gas is air, composed of 
5 species (Ah O, iV 2 , 0o, NO), the free stream Mach 
number is M = 15, the free-stream temperature is 
Toe = 300°/^. Three cases of free-stream pressure 
will be considered, leading to three stiffness condi- 
tions: 

easel : = 10“ 5 atm 

case2 : = 10~ 4 atm 

case3 : P^ = 10~ 2 atm 

The stiffness is defined as the ratio of the largest 
time scale (here presumably the convective one) to 
the smallest (chemistry). The convective time scale 
is obtained from the choice of Courant (CFL) number 
we choose to run the simulation at. The chemistry 
time scale can be defined in two ways: 

• an intrinsic time scale, obtained from the maxi- 
mum rate of change of any chemical specie For 
example, the chemical timescale will be the time 
required for a specie molar fraction to change by 
more than 10%. provided it is not close to zero 
initially. 

• a coupling time scale, defined as the time re- 
quired for the chemistry to modify any flow vari- 
able by (say) more than 5%. Since the chemistry 
affects mostly the temperature, this is the vari- 
able used in that case. 

The second time scale provides a global limitation on 
the time step to be used: if the chemical effects dra- 
matically change the formation energy of the mixture 
during the time step, and if this 6e° is large com- 
pared to the internal energy, the numerical solution 
becomes rapidly unstable. This has a profound effect 
on the choice of numerical methods to be used, for 
example, in combustion. In the remainder of this 
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paper, we always limit the global time step such 
that the estimated relative change in temperature, 
induced by chemical reactions, is smaller or equal to 
5%. We allow the use of large (CFL> 1) global time 
steps provided this condition is satisfied Let us em- 
phasize that this restriction still allows us to consider 
stiff problems, where the stiffness is defined by using 
the intrinsic chemica’ time scale Chemical equilib- 
rium can be reached rapidly, without significantly 
modifying the temperature, the flow conditions sim- 
ply must be such that the equilibrium values are not 
very different from the initial values, or that the en- 
ergies of formation are relatively small compared to 
the internal energy. 

The results presented in this section are obtained 
using the following methods^ : 

1 FC/Block tridiagonal 

2 OS/Block tridiagonal 

3 EC/Block tridiagonal 

4 EC/Scalar tridiagonal 

5 EC/N S — Split tridiagonal 

The chemistry is always computed with a sing 1 * it- 
eration. 

Let us look first at case 1, for P » = 10 -5 atm. 
The profiles of temperature are shown (left scale on 
the plots) as well as the mole fractions of N and O 
atoms (right scale on the plots). Figure 2 shows these 
profiles for a calculation at CFL=2 and 4. Figure 2-a 
shows only the results for the FC/Block, EC/Scalar 
and EC/N, -Split methods: the remaining cases of 
EC/Block and OS/Block would show profiles exactly 
identical to the FC/Block method. The agreement 
between the other methods is also quite good. This 
is also true for the CFL=4 case (Figure 2-b, bottom), 
although to a lesser extent: in this figure, the curves 
for the EC/Block case are omitted, since they are 
identical to the FC/Block results. It appears there- 
fore for this case that the EC and OS methods are 
as accurate as the FC method. The N, -Split method 
shows slight errors in species concentrations near the 
shock, in the region of highest concentration gradi- 
ents, which worsen as the CFL number grows. The 
Scalar method has an overshoot at CFL— 2, and can- 
not be operated at larger C^L numbers. All methods 
fail for larger CFL values. 

Figure 3 shows the same profiles (obtained when 
the shock reaches the same position) for a slightly 

tThe not&tior us«d has been mentioned in the previous 
section, and is summarized in Table 1. 


stiffer problem (case 2). Again, we had perfect 
agreement between OS/Block, EC/Block and 
FC/Block methods, and the EC/Block profiles were 
omitted for clarity. The overshoots in mole fraction 
near th* shock, for the Ng -Split method, are worse 
tor this stiffer case. Again, the Scalar method works 
reasonably well for CFL=2. but fails for larger val- 
ues These calculations were done using the standard 
mtnmod limiter in the convective fluxes, as described 
by Harten [3], with an entropy fix € cs 0.1. When us- 
ing a more compressive flux-limiter, such as the ‘Su- 
perbee’ limiter, the calculation could proceed as well, 
although with very slight oscillations. Reducing the 
entropy parameter to i ^ 0.01 would lead to more 
severe oscillations. Therefore the rapid elimination 
of the transients can be best achieved by ensuring 
that sufficient numerical diffusion is present. The fi- 
nal flow solution therefore would need *o be further 
sharpened, when the steady state is nearly achieved. 

Figure 4 shows the stiffest case for CFL=2. All 
methods failed for larger CFL values. It is remark- 
able that the FC/Block method failed for this case, 
while the EC/Block method gives the best results. 
The OS/Block method (which uses the fractional 
step approach) gives very similar results, and can be 
considered as accurate. Surprisingly, the EC/Scalar 
method is stable, although the species profiles show 
an unphysical kink in the relaxation region. In order 
to better determine which method is more accurate, 
we computed the same case on a larger grid (2000 
points) using the FC/Biock method. By increasing 
the grid density, we achieved a ten-fold reduction 
of the stiffness of the problem. The FC/Block was 
then run successfully, and gave a very short relax- 
ation zone (see Figure 4b). We also attempted to 
better reproduce this relaxation on the coarse (200 
points) grid by either 1) sub-cycling the chemistry 
or 2) reducing the time step. Figure 4b shows the 
comparison, for example between the FC/Block re- 
sults computed on the high-density grid, with the 
EC/Scalar with 10 sub-iterations of the chemistry 
and an explicit calculation (CFL=0.2). The two lat- 
ter cases are not very different from the results of 
Figure 4a, i.e., neither the sub-iterations nor the time 
step reduction greatly improve the solution. It seems 
that all methods tend to overestimate the length of 
the chemical relaxation zone in stiff cases, although 
the final equilibrium result is accurate. We must 
conclude also that the EC or OS methods are more 
stable than the FC method in stiff cases: we also ob- 
served this feature on other stiff cases. The mixing 
of non-diagonal elements in the global Jacobians, be- 
tween convective and chemical terms, may make the 
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matrices more prone to ill-conditionality, »"d reduce 
the stability 

It seems therefore that only relatively small \ a: 
ues of the CFL number can be effectively us;d for t. • 
transients, and therefore only the Scalar me* hod, so 
far, is efficient. However, it is not accurate enough 
when the chemistry is very stiff. It also appears that 
the .V, -Split method, in its current form, suffers from 
unphysical numerical species diffusion in the region 
of strong gradients, and for large time steps (this er- 
ror is inexistent in the explicit regime). Since there 
are other cases where implicit methods can have a 
significant impact, we will look also at expanding 
flows in the next section. 

IV. One-Dimensional Nozzle 

We will model a converging/diverging nozzle, 
with 150 grid points in the axial direction: the calcu- 
lations were performed using two-dimensional codes, 
and the grid used 10 points in the radial direction. 
Since we were interested only into axial profiles, this 
was considered sufficient for our purposes. Notice 
that now there will be an additional error due to 
the dimensional split in the implicit methods. The 
left boundary condition and initial state considered 
a gas at a pressure of 4.205 atm and 1000°K. How- 
ever, the gas composition was arbitrarily set to non 
equilibrium values by increasing the amount of disso- 
ciation: this had the effect of stiffening the chemistry 
in the subsonic region of the nozzle. The calculations 
were always started impulsively, and run at CFL< 1 
(explicit) until the shock exits the nozzle, before the 
implicit models were tried. We used this case to eval- 
uate the effect of sub-iterations and sub-coupling in 
the chemistry. 

Figure 5 shows the comparison of residual his- 
tory for the 4 implicit methods used, i.e. FC/Block, 
EC/Block, EC/Scalar and EC/N, -Split, and with- 
out any sub-iterations or sub-coupling. The resid- 
ual of the subsonic zone (solid line) and supersonic 
zone (dotted line) have been shown separately. A 
first break point in the curves shows the end of the 
explicit pre-calculation, used for elimination of the 
shock from the nozzle. The implicit scheme is then 
used, with a constant CFL=1.5, until a time of 5 mil- 
liseconds. At that point, the calculation is pursued 
further for the supersonic region only, the subsonic 
region remaining frozen. The residual shown is for 
the total energy, and is averaged over the entire vol- 
ume of the region considered. 


Figure 5-a shows the results for the FC/Block 
method, Figure 5-b for the EC/Block method. Both 
show good convergence properties, with a slight im- 
provement for the EC/Biock method. The A r , -Split 
method (Figure 5-c) hae- even better convergence 
properties, but the Scalar method (Figure 5-d) shows 
a non-vanishing residual for the subsonic region All 
methods converge rapidly in the supersonic region 
when computed separately, indicating that most of 
the problems (if any ) are located in the subsonic re- 
gion. 

The solutions obtained at 5 milliseconds are plot- 
ted in Figure 6, for the atomic oxygen mole fraction. 
The solution for the scalar method is slightly in error 
in the subsonic region, but quickly recovers during 
the expansion and leads to the correct final value 
The N, -Split method has the opposite behavior, i.e. 
has an error increasing with the distance along the 
nozzle: it seems that the species convection suffers 
from some unphysical diffusion of species, also no- 
ticed in the results of the previous section: there 
is a phase error between each species convection, 
which is irreversible. By contrast, the Scalar method 
correctly propagates the species, but does not ac- 
curately couple the convection with the momentum 
and energy equations. This may lead to fluctuations 
in pressure or temperature, which quickly disappear 
when the flow becomes near supersonic. 

The use of sub-iterations in the chemistry did 
not change the results for this case. Increasing the 
stagnation pressure and the stiffness slowly lead to 
noticeable effects. The most dramatic differences 
between the cases of sub-iterated and non-iterated 
chemistry were observed for very stiff systems, at 
the limit of stability. In order to demonstrate the 
effect of sub-iterations, or sub-cycling of the chem- 
istry, we consider a high pressure (400 atmospheres) 
case, with an initial temperature of 6000 °K , and 
a highly non-equilibrium initial composition (non- 
dissociated air). A constant time step of 5 x 10~ 6 
seconds was assumed. Figure 7 shows the results 
of the chemistry integration (no fluid dynamics) for 
both non-iterated and sub-iterated (20 cycles) cases. 
It is clear that a single step of the chemical integra- 
tion, with At = 5 x 10~ 8 sec, leads to very large 
changes in spe^es concentrations and temperature. 
This will significantly affect the remainder of the his- 
tory of the chemical integration. If the time step is 
not too large, the correct equilibrium values may be 
obtained in the final steady state: if the time step 
is large enough, unphysical values (i.e. negative con- 
centrations) may be obtained during the first step, 
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and the correct solution cannot be recovered. By 
using a smaller time step (10~ d sec), both methods 
give essentially the same history, and the same fi- 
nal values, which agree with the values obtained in 
Figure 7 for the sub-iterated case It is clear then 
that in some severe cases, the sub-iteration of the 
chemistry can yield a higher stability and a higher 
accuracy: these cases may be found for example in 
high-pressure shocks, detonations, or strong ionizing 
shocks, when the grid used *s coarse. 

The method of sub-coupling was also tested on 
some other stiff cases. It was found that the stabil- 
ity was slightly reduced when the sub-coupling was 
incorporated. When the chemistry is sub-iterated 
and sub-coupled, the convection of species during the 
sub-step St is estimated and included in the varia- 
tion. While this process accounts for the effect of 
the convection on the chemistry, it fails to account 
for the reverse process, and it fails to account for 
the influence of other convective terms, specifically 
the pressure waves. The coupling of the chemistry is 
not ‘in phase’ with all the convective equations. In 
subsonic regions and behind shock waves, the pres- 
sure waves are a dominant j. . ocess, and a signifi- 
cant error is made. We therefore recommend that 
no sub-coupling (SIC) be used if the chemistry is 
sub-iterated. 

V: Two-Dimensional Shock Layer 

The final test will be done for a 2-dimensional, 
axi-symmetric flow around a blunt body. This is a 
typical flow configuration of interest. The flow is 
modeled with a 144x80 grid, the free stream Mach 
number is Mqq — 25, the free stream tempera- 
ture and pressure are T « = 241.75°K and P ^ = 
1.65 10“ 4 atm; the free stream is air, modeled using 
11 species TV, 0 ^ 202 ^ 0 ^ + , 0 + y N+ , 0 + ,N 0 + 
and a 15 reaction set from Dunn ic Kang [5]. The 
blunt body shape is taken from the Apollo space- 
craft. 

The calculations were proceeded with the 
FC/BIock method and the OS/Scalar method. Since 
the chemical charges were quite important in the 
shock layer, our restriction in time step due to the 
temperature changes (relative change < 5%, see sec- 
tion III) prevented us to compute the flow implicitely. 
Any attempt to increase the time step, and therefore 
tc allow a larger change in temperature due to chem- 
ical effects led very quickly to flow instabilities. The 
comparison presented below is therefore between a 
FC/Block implicit method run at small time steps, 


and an Operator-Split method where the fluid dy- 
namics are computed expliciteiy. The pressure be- 
hind the shock is close to 0 5 atmospheres, with a 
temperature between 8,000 and 12,000 ®K, and the 
chemistry is rather stiff, especially due to the reac- 
tions involving electrons Figure 8 shows the compar- 
ison between the two methods along the stagnation 
line, and the agreement is very satisfactory. 

There is a lack of resolution of the shock front, 
and we proceeded to improve the results by adapting 
the grid in the neighborhood of the shock. Several 
adaptions were performed, first on the temperature 
gradients, then on the chemical gradients (;V 2 ). Af- 
ter each adaption, the flow was computed further 
until convergence. The adaption procedure used the 
SAGE code [6] developed at Ames, and affected grid 
points in a direction normal to the blunt surface only. 
Figure 9 shows the comparison between the origi- 
nal, non-adapted polution and the results from the 
final adaption. Since the results from the FC/ lock 
and OS/Scalar methods were found to be agreement, 
and since the latter method is considerably more ef- 
ficient, only the OS/Scalar method was used for the 
adapted cases. We see in Figure 9 that the peak 
temperature has changed significantly (15 %), and 
so has the shock location. Although the flow vari- 
ables relax to the same values in the midst of the 
shock layer, the unresolved relaxation zone may still 
affect some important engineering variables, such as 
the radiative heating at the wall. The radiative emis- 
sion power behind the shock will depend strongly on 
the temperature and species densities, both varying 
rapidly in that region, and being very sensitive to 
the grid resolution. Therefore, ? radiation code was 
used to compute the intensity along the line-of-^ght 
in the stagnation region, and the heat flux at the 
wall: this computation was performed after the flow 
steady state was obtained, i.e. the flow and radiation 
were not coupled. After each adaption, the change 
in radiative heat flux was computed and compared: 
the results are shown in Table 2, for both the opti- 
cally thin case and the optically thick case. In the 
former case, the relative changes are quite large, and 
the values converge slowly. In the thick case, the 
absorption by the core of the shock layer tends to 
damp the perturbations: for that case, we see that 
the heat flux converges more rapidly towards a fi- 
nal value. Since the relative change is small (-1.8%) 
after the Z rd adaption, we considered t^at the res- 
olution was now sufficient. The comparison in ra- 
diative spectrum at the wall between the unadapted 
and final solution is shown in Figure 10. Most of the 
changes occur in the UV region. 
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It is important to remark that the computed radi- 
ation did not include the continuous spectrum, and 
therefore the variations in radiative heating at the 
wall may be under-estimated Additionally, the den- 
sity is sufficiently large in this example that equi- 
librium radiation can be assumed: this considerably 
reduces the uncertainty in the somputation of the 
radiative flux. For lower density and higher veloc.iy 
cases one must include thermal non-equilibrium ef- 
fects. The relaxation zone becomes then even more 
important to resolve accurately. 

It appears from this example that in practice it 
will be very difficult to compute a flow using an im- 
plicit method with a large CFL number, and that 
hydrodynamics-chemistry coupling effects will some- 
times limit the time step to CFL values below 1. Ad- 
ditionally. we may be required in practice to trans- 
form the grid, according to the solution obtained, in 
order to reach the desired accuracy: these adaptions 
need to be performed several times. It would seem 
therefore that a more efficient approach would com- 
bine the flow computation with the grid refinement 
Indeed, there is a technique that can potentially lead 
to more efficient computations: using unstructured 
grids, the computat'onal cells can be subdivided at 
will to give better accuracy in the regions that need 
it. Similarly, the subdivided cells can be regrouped 
in regions of low gradients, in order to keep the total 
number of cells within reasonable limits. Such a tech- 
nique would use a small number of cells to start with, 
and progressively refine them: most of the transients 
would then be computed using a small number of 
cells, leading to a more efficient procedure. 

VI: Conclusions and Recommendations 

We have not used here all the possible variations 
on the implicit schemes, neither have we exhausted 
the methods of coupling the chemistry and the fluid 
dynamics. We have however used techniques which 
are commonly used, and, we hope, demonstrated the 
trends for practical problems. We can draw several 
conclusions from this study: 

1. It is clear that on many problems of interest, 
the calculations cannot proceed with very large 
CFL numbers during the approach to steady 
state. Inevitably, for large numbers of species, 
the Block Tridiagonal methods cease to be ef- 
ficient in that regime Only Scalar Tridiagonal 
methods, or even explicit methods remain effi- 
cient. 


2. It is clear that the Operator-Splitting approach, 
including the Explicit-Coupling between chem- 
istry and hydrodynamics, is at least as accurate 
as the Fully Coupled, and apparently more sta- 
ble for very stiff problems. Sub-iterations of the 
chemistry can further improve the accuracy and 
stability in the most severe cases of stiffness 

3. The A, -Split method, at least in its present for- 
mulation, is too inaccurate for large time steps 
or strong concentration gradien • *. This disa- 
pointing result is not completely understood at 
the moment. It does not affect our conclusions, 
since the method is less performant than the 
Scalar method This results should however be 
investigated further, since it may have implica- 
tions on other systems, such as a two- or three- 
fluid plasma, where the implicit treatment of the 
electron component gas dynamics is mandatory. 

4. Calculations of shock layers on fixed grids may 
not be sufficiently accurate if radiative phenom- 
ena or thermal non-equilibrium effects must be 
considered. In the example shown, several iter- 
ations at grid adaption were necessary. Other 
calculations on similar problems were also per- 
formed, that supported this conclusion. It ap- 
pears then that dynamical grid adapting should 
be performed during the course of tne calcula- 
tion, for higher efficiency. 

Dynamic grid refinement could lead to even 
higher efficiencies if both the distribution and the 
overall number of grid points are allowed to vary. 
This can be done on structured as well as unstruc- 
tured grids. The construction of implicit schemes on 
unstructured grids would be quite complex. How- 
ever, we have concluded that this may not be a re- 
striction for many cases. An explicit algorithm will 
therefore be sufficient, and the technique is reduced 
to a sophisticated book-keeping problem. In addi- 
tion, the use of explicit, Operator-Split techniques 
allows us to take advantage of massively parallel (or 
mixed) computer architectures. This method will be 
investigated in the future. 

We have net mentioned another technique appli- 
cable for Operator-Split methods, when the chem- 
istry is very stiff. The chemistry (or other inter- 
nal process) can be rescaled, or ‘slowed-down’ ar- 
tificially: this may have the effect of increasing the 
relaxation distances. However, we have made pre- 
liminary calculations that seem to indicate that in 
the very severe cases of stiffness, the changes are not 
perceptible. In addition, this procedure can be used 
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during the elimination of the transients, then the 
rescaling is progressively eliminated until a steady 
state with the proper time scale is obtained. If the 
rescaling is not eliminated, a false steady solution 
is obtained. The influence of the numerical proce- 
dure on the steady solution is also a serious ques- 
tion, discussed recently by Lafon k Yee [7]. They 
show that for flows coupled with non-linear source 
terms, the steady state reached may depend on the 
path used to reach it. It is clear therefore that the er- 
rors induced by the numerical procedures can never 
be under-estimated, and that all users of CFD should 
proceed with extreme caution. 
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designation 

Treatment of hydrodynamics 

Chemistry/Convection coupling method 

FC/Block 
EC/Block 
OS/Block 
EC/Scalar 
EC/N, -Split 
EC-Sl/Block 
EC-SIC/Block 

N x N Block Tridiagonal (section II-A) 
N x N Block TYidiagona! (section II- A) 
N x N Block Tridiagonal (section II-A) 
Scalar Tridiagonal (section Il-B) 

N, -split method (section II-C) 

N x A r Block Tridiagonal 
;V x N Block Tridiagonal 

Implicit (or Fully-) coupled 
Explicit (or Loosely-) coupled 
Operator-Split (or Fractional Step) 

Explicit coupling 
Explicit coupling 

Explicit with Sub-It^ra ions (or Sub-Cycling) 
Explicit. Sub-Iterations and Sub-Coupling 


Table 1: Designation of numerical methods and coupling methods used in this study. 


Grid Cycle 

Relative Change 

Relative Change 


in Surface Flux 

in Surface Flux 


Optically Thin Gas 

Optically Thick Gas 


12000-8000 A] 

[1740-1750 A] 

Non-adapted Grid - Adaption 1 

-25.1 % 

-11.0 7c 

Adaption 1 - Adaption 2 

+27.1 % 

+5.2 % 

1 Adaption 2 - Adaption 3 

-3.5 % 

-1.8 % 


Table 2: Axisymmetric blunt body calculations - results of grid adap n study. 
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Required Memory (MW) Normalized CPU Time (CPU} mp H C jt/CPU e ,ipHdt) 



Figure 1-a. Normalized CPU requirement for Block, Scalar and N s -Split methods 
(nonreacting). CPU is obtained as time per iteration per grid point, and normalized 
to CPU requirement for explicit method. 



Number of Species (N s ) 

Figure 1-b. Memory requirement versus number of speoes for same methods 
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• FC/BIock 

EC/Scalar 

- EC/N, -Split 


-25 -20 -15 

Distance from wall (cm) 

Figure 2-a. 1-D Shock (Case 1, P_=10 s atm). Results for CFL=2. 
Results for EC/Block and OS/Block methods were indistinguishable 
from FC/Block method and hence, are not shown. 
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Figure 2-b. 1-D Shock (Case 1, P_=1Q' S atm.). Results for CFL=4. 
The EC/Scalar method failed for this CFL number. 
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Temperature (xlO 1 K) Temperature (xlO 3 K) 
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EC/N,-Split 
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Distance from wall (cm) 

Figure 3-a. 1-D Shock (Case 2, P_=10 4 atm.). Results for CFL=2. 

Results of EC/Block and OS/Block methods are indistinguishable 
from FC/ Block results and hence, are no> shown. 
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OS/Block 
EC/N, -Split 


Distance from wall (cm) 

Figure 3-b. 1-D Shock (Case 2, P_=10' 4 atm.). Results for C.'L=4 
The EC/Scaiar method failed for this case. 
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EC/Block 

- - OS/Block 
EC/Scalar 


-28 -27 -26 -25 -24 -23 

Distance from wall (cm) 

Figure 4-a. 1-D Shock (Case 3, P.=10’ 2 atm.). Results for CFL=2. 
The FC/Block method failed and the EC /N, -Split method gave 
nonphysical, oscillatory solution. 


FC/Block (High Res.) 

& EC/Scalar/Subcycled (Low Res.) 
-A- EC/Block/Explicit (Low Res.) 



-27.5 


-27.0 


Distance from wall (cm) 


f 18 im e ^ ^* tails ° f Case 3 - FC/Block solution is obtained for a higher density 
(x 10) gnd. Other solutions shown are for the standard grid but with subcycled Y 
chemistry or smaller time steps (CFL=0.2) y ea 
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Temperature (K) Mole fraction 


FC/BIock 

— OS/Block 
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Distance along nozzle axis (m) 

Figure 6. Quasi-ID Nozzle * steady state axial distribution of atomic oxygen. 
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Temperature (Subcycles=20) 

Temperature (Subcycles=l) 

Atomic Oxygen (Subcycles=20) 

-■o - Atomic Oxygen (Subcycles=l) 


120 


b 100 




0.0 


T 

0.2 


T~ 

0.4 


0.6 


77 


b 80 


b 60 


b 40 


20 


1.0 


Time (xlO -6 s) 


f*^ C J-^ aSI ; 1D NoZ2l ,V TeSt of subiterated chemistry (EC/Block method). 
Initial state far from equilibrium and global time step (At) is 50 ns. Subcycled 
case (20 cycles) has an effective time step (50 of 25 ns. The temperature is 
recomputed only at the end of the integration over the global time step 50 ns. 


Molar Concentration (moles/m 3 ) 



Mole Fraction (xlO 3 ) Temperature (xlO 3 K) 



Figure 8-a. Axisymmetric Blunt Body - Temperature 
profile along the stagnation line (non-adapted grid) 



Figure 8-b. Axisymmetric Blunt Body - Electron mole fraction 
profile along the stagnation line (non-adapted grid). 
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Mole Fraction of e' (xlO* ) Temperature (xlO 3 K) 


Non-adapted Grid 
Adapted Grid (N 2 Gradient) 


Distance from body surface (xlO * 3 m) 

Figure 9-a. Axisymmetric Blunt Body * Temperature profile along 
the stagnation line for non-adapted and adapted (final pass) grid. 
Note the significant change in shock location and peak temperature. 


Non-adapted Grid 
Adapted Grid (N 2 Gradient) 



Distance from body surface (xlO m) 

Figure 9-b. Axisymmetric Blunt Body - Electron mole fraction profile 
along stagnation line for non-adapted and adapted (final pass) grids. 
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Wavelength (A) 

Figure 10. Axisymmctric Blunt Body - Comparison of radiation spectra 
incident on the wall (stagnation line of sight). The change is • 35% in 
the UV(2500 A) region. 
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COMPUTATION OF NONEQUILIBRIUM VISCOUS FLOWS 
IN ARC-JET WIND TUNNEL NOZZLES 


by 

Tahir Gok^cn f 

Eloret Institute, 3788 Fabian Way, Palo Alto, CA 94303 


Arc-jet hypersonic wind tunnels are the means for studying aerodynamic heating and 
ablation associated with hypersonic flight. In these tunnels a high enthalpy gas produced 
m the arc chamber, vibrationally excited and dissociated, is expanded through a nozzle 
into the test section. Analysis of the nozzle flow is necessary for calibration of the wind 
tunnel flow. Typically these nozzle flows can be characterized as follows: an equilibrium 
flow in the converging section of the nozzle; a nonequilibrium flow near the throat; a 
frozen flow downstream in the diverging section. Therefore, a nonequilibrium formulation 
is required for the analysis of entire flow. Such a formulation needs to be validated against 
experimental data in the equilibrium as well as in the nonequilibrium flow regime. The 
nonequilibrium formulation using Park’s two-temperature thermochemical model has pro- 
duced very satisfactory results in dissociation dominated flows. 1 However, it has not been 
evaluated in the expanding nozzle flows where recombination and vibrational de-excitation 
arc the do m i n ant noncquilibrinm processes. 


Objectives of the proposed paper are to present an analysis of two dimensional or 
axisymmetric arc-heated wind tunnel nozzle flows, and assess the accuracy of the nonequi- 
librium formulation by comparing the computed flow quantities with existing experimental 
data. The comprehensive experimental data set of Folk and Smith 1 taken at a direct con- 
nect arc-heated wind tunnel nozzle of Air Force Flight Dynamics Laboratory is suitable 
for this purpose. 

In the present formulation, two dimensional/ axisymmetric nonequilibrium viscous 
flows of air are considered. The governing Navier-Stokes equations are augmented with 
the equations accounting lor thermochemical nonequilibrium processes. The equation set 
consists of fifteen partial differential equations: eleven mass conservation equations for 
species, two momentum equations for two dimensional flows, and two energy equations. 
The present nonequibbnum gas model for air consists of eleven chemical species, ( N 3 , 0 3 , 
NO, N, O, Nf, Of, NO*, N+, 0 + , e~), and the thermal state of the gas is described 
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by two temperatures: translational-rotational and vibrational-electronic. The chemical re- 
action rates are prescribed by Park’s two-temperature model where the basic dissociation 
rate is assumed to be governed by the geometric average of translational and vibrational 
temperatures. 5 The paper will give a more detailed description of the model. The govern- 
ing equations are discretized using a finite volume approach and the resulting difference 
equations are solved using an implicit Gauss-Seidel relaxation method. The method is 
fully implicit for fluid dynamics and thermochemistry and uses flux vector splitting for 
convective fluxes. Details of the method can be found in Refs. 4-6. 

In the proposed paper, evaluation of the nonequilibrium formulation will be made 
through comparisons of the computed flow quantities at the nozzle exit with the exper- 
imental data of Folk and Smith. 3 For a wide range of reservoir conditions, exit surveys 
and centerline measurements of Pitot pressure, mass flux, stagnation point heat transfer 
rate, and wall static pressure of a 2-ft-diameter conical nozzle are available. Since the wall 
pressure and radial mass flux measurements at the nozzle exit cover the nonequilibrium 
as well as the equilibrium regimes, assessment of nonequilibrium thermochemical models 
will be possible by comparing the computed results with the experiment. Although the 
wall pressure and mass flux measurements can be compared directly with computations, 
comparison of the parameters such as Pitot pressure and stagnation point heat transfer 
requires a correlation between flow quantities and these parameters. This will be achieved 
by computational simulation of the flow at the stagnation point of a sphere (or the probe 
itself) or by using expressions such as the one by Fay and Riddell. 7 The computed total 
enthalpy profiles will be compared with the experiment using the stagnation point heat 
transfer data. These comparisons are needed for the calibration of arc-jet nozzle flows. 

As a representative case, an axi symmetric nonequilibrium flow in a converging- 
diverging conical nozzle is computed. The initial conditions at the nozzle entrance are 
specified as air at temperature of T 0 = T v0 = 5000 °K and pressure of p 0 = 35 atm. The 
exit pressure is specified such that shockless supersonic flow is sustained in the diverging 
section of the nozzle. The nozzle wall is assumed to be at a constant temperature of 
T w = 500 °K and non-catalytic to recombination. The nozzle has conical angles of 30 ° for 
converging and 16 ° for diverging section, and its entrance, throat and exit diameters are 
20, 2.54, and 63.5 cm, respectively. This particular geometry is taken from Ref. 2. Figure 1 
depicts the nozzle geometry and computational grid. For the reservoir conditions chosen, 
the flow remains in thermochemical equilibrium in converging section and goes through 
a nonequilibrium expansion downstream of the nozzle throat. Recombination and vibra- 
tional de-exdtation of air are principal nonequilibrium processes in the nozzle expansion. 
Figures 2 and 3 show the flowfield contours of atomic oxygen and vibrational-electronic 
temperature. In Fig. 2, recombination of oxygen downstream of the throat and in the 
viscous boundary layer is clearly observed. Also it is interesting to note that in Fig. 3, the 
vibrational temperature rise in the diverging section of the nozzle after dropping near the 
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throat. This is mainly due to atomic recombination. In Fig, 4, variation of the computed 
temperatures along the nozzle centerline is presented. Figure 5 shows radial variation of 
computed temperature profiles at the nozzle exit where the temperatures are normalized 
by corresponding centerline values. 

Since the computed results are preliminary, no comparisons with experimental data 
is presented. The proposed paper will make detailed comparisons and assess the accuracy 
of thermochemical nonequilibrium models. 
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